/*
 * Copyright 2012-2017 Jeff Hain
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
/*
 * =============================================================================
 * Notice of fdlibm package this program is partially derived from:
 *
 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
 *
 * Developed at SunSoft, a Sun Microsystems, Inc. business.
 * Permission to use, copy, modify, and distribute this
 * software is freely granted, provided that this notice 
 * is preserved.
 * =============================================================================
 */
package com.nulldev.util.data.FastMath.impls;

import com.nulldev.util.VariableAPI.RandomUtil;

/**
 * Faster (hopefully) versions of java.lang.Math methods, plus additional ones.
 * Cf. README.txt for more info.
 */
@SuppressWarnings("unused")
public final class Jafama extends com.nulldev.util.data.FastMath.Math {

	public static final Jafama INSTANCE = new Jafama();

	/*
	 * For trigonometric functions, use of look-up tables and Taylor-Lagrange
	 * formula with 4 derivatives (more take longer to compute and don't add much
	 * accuracy, less require larger tables (which use more memory, take more time
	 * to initialize, and are slower to access (at least on the machine they were
	 * developed on))).
	 * 
	 * For angles reduction of cos/sin/tan functions: - for small values, instead of
	 * reducing angles, and then computing the best index for look-up tables, we
	 * compute this index right away, and use it for reduction, - for large values,
	 * treatments derived from fdlibm package are used, as done in java.lang.Math.
	 * They are faster but still "slow", so if you work with large numbers and need
	 * speed over accuracy for them, you might want to use normalizeXXXFast
	 * treatments before your function, or modify cos/sin/tan so that they call the
	 * fast normalization treatments instead of the accurate ones. NB: If an angle
	 * is huge (like PI*1e20), in double precision format its last digits are zeros,
	 * which most likely is not the case for the intended value, and doing an
	 * accurate reduction on a very inaccurate value is most likely pointless. But
	 * it gives some sort of coherence that could be needed in some cases.
	 * 
	 * Multiplication on double appears to be about as fast (or not much slower)
	 * than call to <double_array>[<index>], and regrouping some doubles in a
	 * private class, to use index only once, does not seem to speed things up, so:
	 * - for uniformly tabulated values, to retrieve the parameter corresponding to
	 * an index, we recompute it rather than using an array to store it, - for
	 * cos/sin, we recompute derivatives divided by (multiplied by inverse of)
	 * factorial each time, rather than storing them in arrays.
	 * 
	 * Lengths of look-up tables are usually of the form 2^n+1, for their values to
	 * be of the form (<a_constant> * k/2^n, k in 0 .. 2^n), so that particular
	 * values (PI/2, etc.) are "exactly" computed, as well as for other reasons.
	 * 
	 * Tables are put in specific inner classes, to be lazily initialized. Always
	 * doing strict tables initialization, even if StrictFastMath delegates to
	 * StrictMath and doesn't use tables, which makes tables initialization a bit
	 * slower but code simpler. Using redefined pure Java treatments during tables
	 * initialization, instead of Math or StrictMath ones (even asin(double)), can
	 * be very slow, because class loading is likely not to be optimized.
	 * 
	 * Most math treatments I could find on the web, including "fast" ones, usually
	 * take care of special cases (NaN, etc.) at the beginning, and then deal with
	 * the general case, which adds a useless overhead for the general (and common)
	 * case. In this class, special cases are only dealt with when needed, and if
	 * the general case does not already handle them.
	 */

	/*
	 * Regarding strictfp-ness:
	 * 
	 * Switching from/to strictfp has some overhead, so we try to only strictfp-ize
	 * when needed (or when clueless). Compile-time constants are computed in a
	 * FP-strict way, so no need to make this whole class strictfp.
	 */

	// --------------------------------------------------------------------------
	// CONFIGURATION
	// --------------------------------------------------------------------------

	/*
	 * FastMath
	 */

	static final boolean FM_USE_JDK_MATH = getBooleanProperty("jafama.usejdk", false);

	/**
	 * Used for both FastMath.log(double) and FastMath.log10(double).
	 */
	static final boolean FM_USE_REDEFINED_LOG = getBooleanProperty("jafama.fastlog", false);

	static final boolean FM_USE_REDEFINED_SQRT = getBooleanProperty("jafama.fastsqrt", false);

	/**
	 * Set it to true if FastMath.sqrt(double) is slow (more tables, but less calls
	 * to FastMath.sqrt(double)).
	 */
	static final boolean FM_USE_POWTABS_FOR_ASIN = false;

	/*
	 * StrictFastMath
	 */

	static final boolean SFM_USE_JDK_MATH = getBooleanProperty("jafama.strict.usejdk", false);

	/**
	 * Used for both StrictFastMath.log(double) and StrictFastMath.log10(double).
	 * True by default because the StrictMath implementations can be slow.
	 */
	static final boolean SFM_USE_REDEFINED_LOG = getBooleanProperty("jafama.strict.fastlog", true);

	static final boolean SFM_USE_REDEFINED_SQRT = getBooleanProperty("jafama.strict.fastsqrt", false);

	/**
	 * Set it to true if StrictFastMath.sqrt(double) is slow (more tables, but less
	 * calls to StrictFastMath.sqrt(double)).
	 */
	static final boolean SFM_USE_POWTABS_FOR_ASIN = false;

	/*
	 * Common to FastMath and StrictFastMath.
	 */

	/**
	 * Using two pow tab can just make things barely faster, and could relatively
	 * hurt in case of cache-misses, especially for methods that otherwise wouldn't
	 * rely on any tab, so we don't use it.
	 */
	static final boolean USE_TWO_POW_TAB = false;

	/**
	 * Because on some architectures, some casts can be slow, especially for large
	 * values. Might make things a bit slower for latest architectures, but not as
	 * much as it makes them faster for older ones.
	 */
	static final boolean ANTI_SLOW_CASTS = true;

	/**
	 * If some methods get JIT-optimized, they might crash if they contain "(var ==
	 * xxx)" with var being NaN (can happen with Java 6u29).
	 * 
	 * The crash does not happen if we replace "==" with "<" or ">".
	 * 
	 * Only the code that has been observed to trigger the bug has been modified.
	 */
	static final boolean ANTI_JIT_OPTIM_CRASH_ON_NAN = true;

	// --------------------------------------------------------------------------
	// GENERAL CONSTANTS
	// --------------------------------------------------------------------------

	/**
	 * Closest double approximation of e.
	 */
	public static final double E = Math.E;

	/**
	 * Closest double approximation of pi, which is inferior to mathematical pi: pi
	 * ~= 3.14159265358979323846... PI ~= 3.141592653589793
	 */
	public static final double PI = Math.PI;

	/**
	 * High double approximation of pi, which is further from pi than the low
	 * approximation PI: pi ~= 3.14159265358979323846... PI ~= 3.141592653589793
	 * PI_SUP ~= 3.1415926535897936
	 */
	public static final double PI_SUP = Double.longBitsToDouble(Double.doubleToRawLongBits(Math.PI) + 1);

	static final double ONE_DIV_F2 = 1 / 2.0;
	static final double ONE_DIV_F3 = 1 / 6.0;
	static final double ONE_DIV_F4 = 1 / 24.0;

	static final float TWO_POW_23_F = (float) NumbersUtils.twoPow(23);

	static final double TWO_POW_24 = NumbersUtils.twoPow(24);
	private static final double TWO_POW_N24 = NumbersUtils.twoPow(-24);

	static final double TWO_POW_26 = NumbersUtils.twoPow(26);
	static final double TWO_POW_N26 = NumbersUtils.twoPow(-26);

	// First double value (from zero) such as (value+-1/value == value).
	static final double TWO_POW_27 = NumbersUtils.twoPow(27);
	static final double TWO_POW_N27 = NumbersUtils.twoPow(-27);

	static final double TWO_POW_N28 = NumbersUtils.twoPow(-28);

	static final double TWO_POW_52 = NumbersUtils.twoPow(52);

	static final double TWO_POW_N55 = NumbersUtils.twoPow(-55);

	static final double TWO_POW_66 = NumbersUtils.twoPow(66);

	static final double TWO_POW_512 = NumbersUtils.twoPow(512);
	static final double TWO_POW_N512 = NumbersUtils.twoPow(-512);

	/**
	 * Double.MIN_NORMAL since Java 6.
	 */
	static final double DOUBLE_MIN_NORMAL = Double.longBitsToDouble(0x0010000000000000L); // 2.2250738585072014E-308

	// Not storing float/double mantissa size in constants,
	// for 23 and 52 are shorter to read and more
	// bitwise-explicit than some constant's name.

	static final int MIN_DOUBLE_EXPONENT = -1074;
	static final int MIN_DOUBLE_NORMAL_EXPONENT = -1022;
	static final int MAX_DOUBLE_EXPONENT = 1023;

	static final int MIN_FLOAT_NORMAL_EXPONENT = -126;
	static final int MAX_FLOAT_EXPONENT = 127;

	private static final double SQRT_2 = StrictMath.sqrt(2.0);

	static final double LOG_2 = StrictMath.log(2.0);
	static final double LOG_TWO_POW_27 = StrictMath.log(TWO_POW_27);
	static final double LOG_DOUBLE_MAX_VALUE = StrictMath.log(Double.MAX_VALUE);

	static final double INV_LOG_10 = 1.0 / StrictMath.log(10.0);

	static final double DOUBLE_BEFORE_60 = Double.longBitsToDouble(Double.doubleToRawLongBits(60.0) - 1);

	// --------------------------------------------------------------------------
	// CONSTANTS FOR NORMALIZATIONS
	// --------------------------------------------------------------------------

	/**
	 * Table of constants for 1/(PI/2), 282 Hex digits (enough for normalizing
	 * doubles). 1/(PI/2) approximation = sum of TWO_OVER_PI_TAB[i]*2^(-24*(i+1)).
	 * 
	 * double and not int, to avoid int-to-double cast during computations.
	 */
	private static final double TWO_OVER_PI_TAB[] =
		{ 0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xe00649, 0x2EEA09,
				0xD1921C, 0xFE1DEB, 0x1CB129, 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 0x3991d6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928,
				0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF, 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 0xF17B3D,
				0x0739F7, 0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 0x91615E, 0xE61B08, 0x659985,
				0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B };

	/*
	 * Constants for PI/2. Only the 23 most significant bits of each mantissa are
	 * used. 2*PI approximation = sum of TWOPI_TAB<i>.
	 */
	private static final double PIO2_TAB0 = Double.longBitsToDouble(0x3FF921FB40000000L);
	private static final double PIO2_TAB1 = Double.longBitsToDouble(0x3E74442D00000000L);
	private static final double PIO2_TAB2 = Double.longBitsToDouble(0x3CF8469880000000L);
	private static final double PIO2_TAB3 = Double.longBitsToDouble(0x3B78CC5160000000L);
	private static final double PIO2_TAB4 = Double.longBitsToDouble(0x39F01B8380000000L);
	private static final double PIO2_TAB5 = Double.longBitsToDouble(0x387A252040000000L);

	static final double PIO2_INV = Double.longBitsToDouble(0x3FE45F306DC9C883L); // 6.36619772367581382433e-01 53 bits of 2/pi
	static final double PIO2_HI = Double.longBitsToDouble(0x3FF921FB54400000L); // 1.57079632673412561417e+00 first 33 bits of pi/2
	static final double PIO2_LO = Double.longBitsToDouble(0x3DD0B4611A626331L); // 6.07710050650619224932e-11 pi/2 - PIO2_HI
	static final double PI_INV = PIO2_INV / 2;
	static final double PI_HI = 2 * PIO2_HI;
	static final double PI_LO = 2 * PIO2_LO;
	static final double TWOPI_INV = PIO2_INV / 4;
	static final double TWOPI_HI = 4 * PIO2_HI;
	static final double TWOPI_LO = 4 * PIO2_LO;

	/**
	 * Bit = 0 where quadrant is encoded in remainder bits.
	 */
	private static final long QUADRANT_BITS_0_MASK = 0xCFFFFFFFFFFFFFFFL;

	/**
	 * Remainder bits where quadrant is encoded, 0 elsewhere.
	 */
	private static final long QUADRANT_PLACE_BITS = 0x3000000000000000L;

	/**
	 * fdlibm uses 2^19*PI/2 here. With 2^18*PI/2 we would be more accurate, for
	 * example when normalizing 822245.903631403, which is close to 2^19*PI/2, but
	 * we are still in our accuracy tolerance with fdlibm's value (but not
	 * 2^20*PI/2) so we stick to it, to help being faster than (Strict)Math for
	 * values in [2^18*PI/2,2^19*PI/2].
	 * 
	 * For tests, can use a smaller value, for heavy remainder not to only be used
	 * with huge values.
	 */
	static final double NORMALIZE_ANGLE_MAX_MEDIUM_DOUBLE_PIO2 = StrictMath.pow(2.0, 19.0) * (Math.PI / 2);

	/**
	 * 2*Math.PI, normalized into [-PI,PI], as returned by
	 * StrictMath.asin(StrictMath.sin(2*Math.PI)) (asin behaves as identity for
	 * this).
	 * 
	 * NB: NumbersUtils.minus2PI(2*Math.PI) returns -2.449293598153844E-16, which is
	 * different due to not using an accurate enough definition of PI.
	 */
	static final double TWO_MATH_PI_IN_MINUS_PI_PI = -2.4492935982947064E-16;

	// --------------------------------------------------------------------------
	// CONSTANTS AND TABLES FOR SIN AND COS
	// --------------------------------------------------------------------------

	static final int SIN_COS_TABS_SIZE = (1 << getTabSizePower(11)) + 1;
	static final double SIN_COS_DELTA_HI = TWOPI_HI / (SIN_COS_TABS_SIZE - 1);
	static final double SIN_COS_DELTA_LO = TWOPI_LO / (SIN_COS_TABS_SIZE - 1);
	static final double SIN_COS_INDEXER = 1 / (SIN_COS_DELTA_HI + SIN_COS_DELTA_LO);

	static final class MyTSinCos {
		static final double[] sinTab = new double[SIN_COS_TABS_SIZE];
		static final double[] cosTab = new double[SIN_COS_TABS_SIZE];
		static {
			init();
		}

		private static strictfp void init() {
			final int SIN_COS_PI_INDEX = (SIN_COS_TABS_SIZE - 1) / 2;
			final int SIN_COS_PI_MUL_2_INDEX = 2 * SIN_COS_PI_INDEX;
			final int SIN_COS_PI_MUL_0_5_INDEX = SIN_COS_PI_INDEX / 2;
			final int SIN_COS_PI_MUL_1_5_INDEX = 3 * SIN_COS_PI_INDEX / 2;
			for (int i = 0; i < SIN_COS_TABS_SIZE; i++) {
				// angle: in [0,2*PI] (doesn't seem to help to have it in [-PI,PI]).
				double angle = i * SIN_COS_DELTA_HI + i * SIN_COS_DELTA_LO;
				double sinAngle = StrictMath.sin(angle);
				double cosAngle = StrictMath.cos(angle);
				// For indexes corresponding to zero cosine or sine, we make sure
				// the value is zero and not an epsilon, since each value
				// corresponds to sin-or-cos(i*PI/n), where PI is a more accurate
				// definition of PI than Math.PI.
				// This allows for a much better accuracy for results close to zero.
				if (i == SIN_COS_PI_INDEX) {
					sinAngle = 0.0;
				} else if (i == SIN_COS_PI_MUL_2_INDEX) {
					sinAngle = 0.0;
				} else if (i == SIN_COS_PI_MUL_0_5_INDEX) {
					cosAngle = 0.0;
				} else if (i == SIN_COS_PI_MUL_1_5_INDEX) {
					cosAngle = 0.0;
				}
				sinTab[i] = sinAngle;
				cosTab[i] = cosAngle;
			}
		}
	}

	/**
	 * Max abs value for index-based reduction, above which we use regular angle
	 * normalization. This value must be < (Integer.MAX_VALUE / SIN_COS_INDEXER), to
	 * stay in range of int type. If too high, error gets larger because index-based
	 * reduction doesn't use an accurate enough definition of PI. If too low, and if
	 * we would be using remainder into [-PI,PI] instead of into [-PI/4,PI/4], error
	 * would get larger as well, because remainder would just provide a double,
	 * while index-based reduction is more accurate, using delta from index values
	 * and HI/LO values.
	 */
	static final double SIN_COS_MAX_VALUE_FOR_INT_MODULO = ((Integer.MAX_VALUE >> 9) / SIN_COS_INDEXER) * 0.99;

	// --------------------------------------------------------------------------
	// CONSTANTS AND TABLES FOR TAN
	// --------------------------------------------------------------------------

	// We use the following formula:
	// 1) tan(-x) = -tan(x)
	// 2) tan(x) = 1/tan(PI/2-x)
	// ---> we only have to compute tan(x) on [0,A] with PI/4<=A<PI/2.

	/**
	 * We use indexing past look-up tables, so that indexing information allows for
	 * fast recomputation of angle in [0,PI/2] range.
	 */
	static final int TAN_VIRTUAL_TABS_SIZE = (1 << getTabSizePower(12)) + 1;

	/**
	 * Must be >= 45deg, and supposed to be >= 51.4deg, as fdlibm code is not
	 * supposed to work with values inferior to that (51.4deg is about
	 * (PI/2-Double.longBitsToDouble(0x3FE5942800000000L))).
	 */
	static final double TAN_MAX_VALUE_FOR_TABS = StrictMath.toRadians(77.0);

	static final int TAN_TABS_SIZE = (int) ((TAN_MAX_VALUE_FOR_TABS / (Math.PI / 2)) * (TAN_VIRTUAL_TABS_SIZE - 1)) + 1;
	static final double TAN_DELTA_HI = PIO2_HI / (TAN_VIRTUAL_TABS_SIZE - 1);
	static final double TAN_DELTA_LO = PIO2_LO / (TAN_VIRTUAL_TABS_SIZE - 1);
	static final double TAN_INDEXER = 1 / (TAN_DELTA_HI + TAN_DELTA_LO);

	static final class MyTTan {
		static final double[] tanTab = new double[TAN_TABS_SIZE];
		static final double[] tanDer1DivF1Tab = new double[TAN_TABS_SIZE];
		static final double[] tanDer2DivF2Tab = new double[TAN_TABS_SIZE];
		static final double[] tanDer3DivF3Tab = new double[TAN_TABS_SIZE];
		static final double[] tanDer4DivF4Tab = new double[TAN_TABS_SIZE];
		static {
			init();
		}

		private static strictfp void init() {
			for (int i = 0; i < TAN_TABS_SIZE; i++) {
				// angle: in [0,TAN_MAX_VALUE_FOR_TABS].
				double angle = i * TAN_DELTA_HI + i * TAN_DELTA_LO;
				double sinAngle = StrictMath.sin(angle);
				double cosAngle = StrictMath.cos(angle);
				double cosAngleInv = 1 / cosAngle;
				double cosAngleInv2 = cosAngleInv * cosAngleInv;
				double cosAngleInv3 = cosAngleInv2 * cosAngleInv;
				double cosAngleInv4 = cosAngleInv2 * cosAngleInv2;
				double cosAngleInv5 = cosAngleInv3 * cosAngleInv2;
				tanTab[i] = sinAngle * cosAngleInv;
				tanDer1DivF1Tab[i] = cosAngleInv2;
				tanDer2DivF2Tab[i] = ((2 * sinAngle) * cosAngleInv3) * ONE_DIV_F2;
				tanDer3DivF3Tab[i] = ((2 * (1 + 2 * sinAngle * sinAngle)) * cosAngleInv4) * ONE_DIV_F3;
				tanDer4DivF4Tab[i] = ((8 * sinAngle * (2 + sinAngle * sinAngle)) * cosAngleInv5) * ONE_DIV_F4;
			}
		}
	}

	/**
	 * Max abs value for fast modulo, above which we use regular angle
	 * normalization. This value must be < (Integer.MAX_VALUE / TAN_INDEXER), to
	 * stay in range of int type. If too high, error gets larger because index-based
	 * reduction doesn't use an accurate enough definition of PI. If too low, error
	 * gets larger as well, because we use remainder into [-PI/2,PI/2], just
	 * provides a double, while index-based reduction is more accurate, using delta
	 * from index values and HI/LO values.
	 */
	static final double TAN_MAX_VALUE_FOR_INT_MODULO = (((Integer.MAX_VALUE >> 9) / TAN_INDEXER) * 0.99);

	// --------------------------------------------------------------------------
	// CONSTANTS AND TABLES FOR ACOS, ASIN
	// --------------------------------------------------------------------------

	// We use the following formula:
	// 1) acos(x) = PI/2 - asin(x)
	// 2) asin(-x) = -asin(x)
	// ---> we only have to compute asin(x) on [0,1].
	// For values not close to +-1, we use look-up tables;
	// for values near +-1, we use code derived from fdlibm.

	/**
	 * Supposed to be >= sin(77.2deg), as fdlibm code is supposed to work with
	 * values > 0.975, but seems to work well enough as long as value >= sin(25deg).
	 */
	static final double ASIN_MAX_VALUE_FOR_TABS = StrictMath.sin(StrictMath.toRadians(73.0));

	static final int ASIN_TABS_SIZE = (1 << getTabSizePower(13)) + 1;
	static final double ASIN_DELTA = ASIN_MAX_VALUE_FOR_TABS / (ASIN_TABS_SIZE - 1);
	static final double ASIN_INDEXER = 1 / ASIN_DELTA;

	static final class MyTAsin {
		static final double[] asinTab = new double[ASIN_TABS_SIZE];
		static final double[] asinDer1DivF1Tab = new double[ASIN_TABS_SIZE];
		static final double[] asinDer2DivF2Tab = new double[ASIN_TABS_SIZE];
		static final double[] asinDer3DivF3Tab = new double[ASIN_TABS_SIZE];
		static final double[] asinDer4DivF4Tab = new double[ASIN_TABS_SIZE];
		static {
			init();
		}

		private static strictfp void init() {
			for (int i = 0; i < ASIN_TABS_SIZE; i++) {
				// x: in [0,ASIN_MAX_VALUE_FOR_TABS].
				double x = i * ASIN_DELTA;
				double oneMinusXSqInv = 1 / (1 - x * x);
				double oneMinusXSqInv0_5 = StrictMath.sqrt(oneMinusXSqInv);
				double oneMinusXSqInv1_5 = oneMinusXSqInv0_5 * oneMinusXSqInv;
				double oneMinusXSqInv2_5 = oneMinusXSqInv1_5 * oneMinusXSqInv;
				double oneMinusXSqInv3_5 = oneMinusXSqInv2_5 * oneMinusXSqInv;
				asinTab[i] = StrictMath.asin(x);
				asinDer1DivF1Tab[i] = oneMinusXSqInv0_5;
				asinDer2DivF2Tab[i] = (x * oneMinusXSqInv1_5) * ONE_DIV_F2;
				asinDer3DivF3Tab[i] = ((1 + 2 * x * x) * oneMinusXSqInv2_5) * ONE_DIV_F3;
				asinDer4DivF4Tab[i] = ((5 + 2 * x * (2 + x * (5 - 2 * x))) * oneMinusXSqInv3_5) * ONE_DIV_F4;
			}
		}
	}

	static final double ASIN_MAX_VALUE_FOR_POWTABS = StrictMath.sin(StrictMath.toRadians(88.6));
	static final int ASIN_POWTABS_POWER = 84;

	static final double ASIN_POWTABS_ONE_DIV_MAX_VALUE = 1 / ASIN_MAX_VALUE_FOR_POWTABS;
	static final int ASIN_POWTABS_SIZE = (FM_USE_POWTABS_FOR_ASIN || SFM_USE_POWTABS_FOR_ASIN) ? (1 << getTabSizePower(12)) + 1 : 0;
	static final int ASIN_POWTABS_SIZE_MINUS_ONE = ASIN_POWTABS_SIZE - 1;

	static final class MyTAsinPow {
		static final double[] asinParamPowTab = new double[ASIN_POWTABS_SIZE];
		static final double[] asinPowTab = new double[ASIN_POWTABS_SIZE];
		static final double[] asinDer1DivF1PowTab = new double[ASIN_POWTABS_SIZE];
		static final double[] asinDer2DivF2PowTab = new double[ASIN_POWTABS_SIZE];
		static final double[] asinDer3DivF3PowTab = new double[ASIN_POWTABS_SIZE];
		static final double[] asinDer4DivF4PowTab = new double[ASIN_POWTABS_SIZE];
		static {
			init();
		}

		private static strictfp void init() {
			if (FM_USE_POWTABS_FOR_ASIN || SFM_USE_POWTABS_FOR_ASIN) {
				for (int i = 0; i < ASIN_POWTABS_SIZE; i++) {
					// x: in [0,ASIN_MAX_VALUE_FOR_POWTABS].
					double x = StrictMath.pow(i * (1.0 / ASIN_POWTABS_SIZE_MINUS_ONE), 1.0 / ASIN_POWTABS_POWER) * ASIN_MAX_VALUE_FOR_POWTABS;
					double oneMinusXSqInv = 1 / (1 - x * x);
					double oneMinusXSqInv0_5 = StrictMath.sqrt(oneMinusXSqInv);
					double oneMinusXSqInv1_5 = oneMinusXSqInv0_5 * oneMinusXSqInv;
					double oneMinusXSqInv2_5 = oneMinusXSqInv1_5 * oneMinusXSqInv;
					double oneMinusXSqInv3_5 = oneMinusXSqInv2_5 * oneMinusXSqInv;
					asinParamPowTab[i] = x;
					asinPowTab[i] = StrictMath.asin(x);
					asinDer1DivF1PowTab[i] = oneMinusXSqInv0_5;
					asinDer2DivF2PowTab[i] = (x * oneMinusXSqInv1_5) * ONE_DIV_F2;
					asinDer3DivF3PowTab[i] = ((1 + 2 * x * x) * oneMinusXSqInv2_5) * ONE_DIV_F3;
					asinDer4DivF4PowTab[i] = ((5 + 2 * x * (2 + x * (5 - 2 * x))) * oneMinusXSqInv3_5) * ONE_DIV_F4;
				}
			}
		}
	}

	static final double ASIN_PIO2_HI = Double.longBitsToDouble(0x3FF921FB54442D18L); // 1.57079632679489655800e+00
	static final double ASIN_PIO2_LO = Double.longBitsToDouble(0x3C91A62633145C07L); // 6.12323399573676603587e-17
	static final double ASIN_PS0 = Double.longBitsToDouble(0x3fc5555555555555L); // 1.66666666666666657415e-01
	static final double ASIN_PS1 = Double.longBitsToDouble(0xbfd4d61203eb6f7dL); // -3.25565818622400915405e-01
	static final double ASIN_PS2 = Double.longBitsToDouble(0x3fc9c1550e884455L); // 2.01212532134862925881e-01
	static final double ASIN_PS3 = Double.longBitsToDouble(0xbfa48228b5688f3bL); // -4.00555345006794114027e-02
	static final double ASIN_PS4 = Double.longBitsToDouble(0x3f49efe07501b288L); // 7.91534994289814532176e-04
	static final double ASIN_PS5 = Double.longBitsToDouble(0x3f023de10dfdf709L); // 3.47933107596021167570e-05
	static final double ASIN_QS1 = Double.longBitsToDouble(0xc0033a271c8a2d4bL); // -2.40339491173441421878e+00
	static final double ASIN_QS2 = Double.longBitsToDouble(0x40002ae59c598ac8L); // 2.02094576023350569471e+00
	static final double ASIN_QS3 = Double.longBitsToDouble(0xbfe6066c1b8d0159L); // -6.88283971605453293030e-01
	static final double ASIN_QS4 = Double.longBitsToDouble(0x3fb3b8c5b12e9282L); // 7.70381505559019352791e-02

	// --------------------------------------------------------------------------
	// CONSTANTS AND TABLES FOR ATAN
	// --------------------------------------------------------------------------

	// We use the formula atan(-x) = -atan(x)
	// ---> we only have to compute atan(x) on [0,+Infinity[.
	// For values corresponding to angles not close to +-PI/2, we use look-up
	// tables;
	// for values corresponding to angles near +-PI/2, we use code derived from
	// fdlibm.

	/**
	 * Supposed to be >= tan(67.7deg), as fdlibm code is supposed to work with
	 * values > 2.4375.
	 */
	static final double ATAN_MAX_VALUE_FOR_TABS = StrictMath.tan(StrictMath.toRadians(74.0));

	static final int ATAN_TABS_SIZE = (1 << getTabSizePower(12)) + 1;
	static final double ATAN_DELTA = ATAN_MAX_VALUE_FOR_TABS / (ATAN_TABS_SIZE - 1);
	static final double ATAN_INDEXER = 1 / ATAN_DELTA;

	static final class MyTAtan {
		static final double[] atanTab = new double[ATAN_TABS_SIZE];
		static final double[] atanDer1DivF1Tab = new double[ATAN_TABS_SIZE];
		static final double[] atanDer2DivF2Tab = new double[ATAN_TABS_SIZE];
		static final double[] atanDer3DivF3Tab = new double[ATAN_TABS_SIZE];
		static final double[] atanDer4DivF4Tab = new double[ATAN_TABS_SIZE];
		static {
			init();
		}

		private static strictfp void init() {
			for (int i = 0; i < ATAN_TABS_SIZE; i++) {
				// x: in [0,ATAN_MAX_VALUE_FOR_TABS].
				double x = i * ATAN_DELTA;
				double onePlusXSqInv = 1 / (1 + x * x);
				double onePlusXSqInv2 = onePlusXSqInv * onePlusXSqInv;
				double onePlusXSqInv3 = onePlusXSqInv2 * onePlusXSqInv;
				double onePlusXSqInv4 = onePlusXSqInv2 * onePlusXSqInv2;
				atanTab[i] = StrictMath.atan(x);
				atanDer1DivF1Tab[i] = onePlusXSqInv;
				atanDer2DivF2Tab[i] = (-2 * x * onePlusXSqInv2) * ONE_DIV_F2;
				atanDer3DivF3Tab[i] = ((-2 + 6 * x * x) * onePlusXSqInv3) * ONE_DIV_F3;
				atanDer4DivF4Tab[i] = ((24 * x * (1 - x * x)) * onePlusXSqInv4) * ONE_DIV_F4;
			}
		}
	}

	static final double ATAN_HI3 = Double.longBitsToDouble(0x3ff921fb54442d18L); // 1.57079632679489655800e+00 atan(inf)hi
	static final double ATAN_LO3 = Double.longBitsToDouble(0x3c91a62633145c07L); // 6.12323399573676603587e-17 atan(inf)lo
	static final double ATAN_AT0 = Double.longBitsToDouble(0x3fd555555555550dL); // 3.33333333333329318027e-01
	static final double ATAN_AT1 = Double.longBitsToDouble(0xbfc999999998ebc4L); // -1.99999999998764832476e-01
	static final double ATAN_AT2 = Double.longBitsToDouble(0x3fc24924920083ffL); // 1.42857142725034663711e-01
	static final double ATAN_AT3 = Double.longBitsToDouble(0xbfbc71c6fe231671L); // -1.11111104054623557880e-01
	static final double ATAN_AT4 = Double.longBitsToDouble(0x3fb745cdc54c206eL); // 9.09088713343650656196e-02
	static final double ATAN_AT5 = Double.longBitsToDouble(0xbfb3b0f2af749a6dL); // -7.69187620504482999495e-02
	static final double ATAN_AT6 = Double.longBitsToDouble(0x3fb10d66a0d03d51L); // 6.66107313738753120669e-02
	static final double ATAN_AT7 = Double.longBitsToDouble(0xbfadde2d52defd9aL); // -5.83357013379057348645e-02
	static final double ATAN_AT8 = Double.longBitsToDouble(0x3fa97b4b24760debL); // 4.97687799461593236017e-02
	static final double ATAN_AT9 = Double.longBitsToDouble(0xbfa2b4442c6a6c2fL); // -3.65315727442169155270e-02
	static final double ATAN_AT10 = Double.longBitsToDouble(0x3f90ad3ae322da11L); // 1.62858201153657823623e-02

	// --------------------------------------------------------------------------
	// CONSTANTS AND TABLES FOR TANH
	// --------------------------------------------------------------------------

	/**
	 * Constant found experimentally: StrictMath.tanh(TANH_1_THRESHOLD) = 1,
	 * StrictMath.tanh(nextDown(TANH_1_THRESHOLD)) =
	 * FastMath.tanh(nextDown(TANH_1_THRESHOLD)) < 1.
	 */
	static final double TANH_1_THRESHOLD = 19.061547465398498;

	// --------------------------------------------------------------------------
	// CONSTANTS AND TABLES FOR ASINH AND ACOSH
	// --------------------------------------------------------------------------

	static final double ASINH_LOG1P_THRESHOLD = 0.04;

	/**
	 * sqrt(x*x+-1) should yield higher threshold, but it's enough due to subsequent
	 * log.
	 */
	static final double ASINH_ACOSH_SQRT_ELISION_THRESHOLD = (1 << 24);

	// --------------------------------------------------------------------------
	// CONSTANTS AND TABLES FOR EXP AND EXPM1
	// --------------------------------------------------------------------------

	static final double EXP_OVERFLOW_LIMIT = Double.longBitsToDouble(0x40862E42FEFA39EFL); // 7.09782712893383973096e+02
	static final double EXP_UNDERFLOW_LIMIT = Double.longBitsToDouble(0xC0874910D52D3051L); // -7.45133219101941108420e+02
	static final int EXP_LO_DISTANCE_TO_ZERO_POT = 0;
	static final int EXP_LO_DISTANCE_TO_ZERO = (1 << EXP_LO_DISTANCE_TO_ZERO_POT);
	static final int EXP_LO_TAB_SIZE_POT = getTabSizePower(11);
	static final int EXP_LO_TAB_SIZE = (1 << EXP_LO_TAB_SIZE_POT) + 1;
	static final int EXP_LO_TAB_MID_INDEX = ((EXP_LO_TAB_SIZE - 1) / 2);
	static final int EXP_LO_INDEXING = EXP_LO_TAB_MID_INDEX / EXP_LO_DISTANCE_TO_ZERO;
	static final int EXP_LO_INDEXING_DIV_SHIFT = EXP_LO_TAB_SIZE_POT - 1 - EXP_LO_DISTANCE_TO_ZERO_POT;

	static final class MyTExp {
		static final double[] expHiTab = new double[1 + (int) EXP_OVERFLOW_LIMIT - (int) EXP_UNDERFLOW_LIMIT];
		static final double[] expLoPosTab = new double[EXP_LO_TAB_SIZE];
		static final double[] expLoNegTab = new double[EXP_LO_TAB_SIZE];
		static {
			init();
		}

		private static strictfp void init() {
			for (int i = (int) EXP_UNDERFLOW_LIMIT; i <= (int) EXP_OVERFLOW_LIMIT; i++) {
				expHiTab[i - (int) EXP_UNDERFLOW_LIMIT] = StrictMath.exp(i);
			}
			for (int i = 0; i < EXP_LO_TAB_SIZE; i++) {
				// x: in [-EXPM1_DISTANCE_TO_ZERO,EXPM1_DISTANCE_TO_ZERO].
				double x = -EXP_LO_DISTANCE_TO_ZERO + i / (double) EXP_LO_INDEXING;
				// exp(x)
				expLoPosTab[i] = StrictMath.exp(x);
				// 1-exp(-x), accurately computed
				expLoNegTab[i] = -StrictMath.expm1(-x);
			}
		}
	}

	// --------------------------------------------------------------------------
	// CONSTANTS AND TABLES FOR LOG AND LOG1P
	// --------------------------------------------------------------------------

	static final int LOG_BITS = getTabSizePower(12);
	static final int LOG_TAB_SIZE = (1 << LOG_BITS);

	static final class MyTLog {
		static final double[] logXLogTab = new double[LOG_TAB_SIZE];
		static final double[] logXTab = new double[LOG_TAB_SIZE];
		static final double[] logXInvTab = new double[LOG_TAB_SIZE];
		static {
			init();
		}

		private static strictfp void init() {
			for (int i = 0; i < LOG_TAB_SIZE; i++) {
				// Exact to use inverse of tab size, since it is a power of two.
				double x = 1 + i * (1.0 / LOG_TAB_SIZE);
				logXLogTab[i] = StrictMath.log(x);
				logXTab[i] = x;
				logXInvTab[i] = 1 / x;
			}
		}
	}

	// --------------------------------------------------------------------------
	// TABLE FOR POWERS OF TWO
	// --------------------------------------------------------------------------

	static final int TWO_POW_TAB_SIZE = USE_TWO_POW_TAB ? (MAX_DOUBLE_EXPONENT - MIN_DOUBLE_EXPONENT) + 1 : 0;

	static final class MyTTwoPow {
		static final double[] twoPowTab = new double[TWO_POW_TAB_SIZE];
		static {
			init();
		}

		private static strictfp void init() {
			if (USE_TWO_POW_TAB) {
				for (int i = MIN_DOUBLE_EXPONENT; i <= MAX_DOUBLE_EXPONENT; i++) {
					twoPowTab[i - MIN_DOUBLE_EXPONENT] = NumbersUtils.twoPow(i);
				}
			}
		}
	}

	// --------------------------------------------------------------------------
	// CONSTANTS AND TABLES FOR SQRT
	// --------------------------------------------------------------------------

	static final int SQRT_LO_BITS = getTabSizePower(12);
	static final int SQRT_LO_TAB_SIZE = (1 << SQRT_LO_BITS);

	static final class MyTSqrt {
		static final double[] sqrtXSqrtHiTab = new double[MAX_DOUBLE_EXPONENT - MIN_DOUBLE_EXPONENT + 1];
		static final double[] sqrtXSqrtLoTab = new double[SQRT_LO_TAB_SIZE];
		static final double[] sqrtSlopeHiTab = new double[MAX_DOUBLE_EXPONENT - MIN_DOUBLE_EXPONENT + 1];
		static final double[] sqrtSlopeLoTab = new double[SQRT_LO_TAB_SIZE];
		static {
			init();
		}

		private static strictfp void init() {
			for (int i = MIN_DOUBLE_EXPONENT; i <= MAX_DOUBLE_EXPONENT; i++) {
				double twoPowExpDiv2 = StrictMath.pow(2.0, i * 0.5);
				sqrtXSqrtHiTab[i - MIN_DOUBLE_EXPONENT] = twoPowExpDiv2 * 0.5; // Half sqrt, to avoid overflows.
				sqrtSlopeHiTab[i - MIN_DOUBLE_EXPONENT] = 1 / twoPowExpDiv2;
			}
			sqrtXSqrtLoTab[0] = 1.0;
			sqrtSlopeLoTab[0] = 1.0;
			final long SQRT_LO_MASK = (0x3FF0000000000000L | (0x000FFFFFFFFFFFFFL >> SQRT_LO_BITS));
			for (int i = 1; i < SQRT_LO_TAB_SIZE; i++) {
				long xBits = SQRT_LO_MASK | (((long) (i - 1)) << (52 - SQRT_LO_BITS));
				double sqrtX = StrictMath.sqrt(Double.longBitsToDouble(xBits));
				sqrtXSqrtLoTab[i] = sqrtX;
				sqrtSlopeLoTab[i] = 1 / sqrtX;
			}
		}
	}

	// --------------------------------------------------------------------------
	// CONSTANTS AND TABLES FOR CBRT
	// --------------------------------------------------------------------------

	static final int CBRT_LO_BITS = getTabSizePower(12);
	static final int CBRT_LO_TAB_SIZE = (1 << CBRT_LO_BITS);

	// For CBRT_LO_BITS = 12:
	// cbrtXCbrtLoTab[0] = 1.0.
	// cbrtXCbrtLoTab[1] = cbrt(1. 000000000000
	// 1111111111111111111111111111111111111111b)
	// cbrtXCbrtLoTab[2] = cbrt(1. 000000000001
	// 1111111111111111111111111111111111111111b)
	// cbrtXCbrtLoTab[3] = cbrt(1. 000000000010
	// 1111111111111111111111111111111111111111b)
	// cbrtXCbrtLoTab[4] = cbrt(1. 000000000011
	// 1111111111111111111111111111111111111111b)
	// etc.
	static final class MyTCbrt {
		static final double[] cbrtXCbrtHiTab = new double[MAX_DOUBLE_EXPONENT - MIN_DOUBLE_EXPONENT + 1];
		static final double[] cbrtXCbrtLoTab = new double[CBRT_LO_TAB_SIZE];
		static final double[] cbrtSlopeHiTab = new double[MAX_DOUBLE_EXPONENT - MIN_DOUBLE_EXPONENT + 1];
		static final double[] cbrtSlopeLoTab = new double[CBRT_LO_TAB_SIZE];
		static {
			init();
		}

		private static strictfp void init() {
			for (int i = MIN_DOUBLE_EXPONENT; i <= MAX_DOUBLE_EXPONENT; i++) {
				double twoPowExpDiv3 = StrictMath.pow(2.0, i * (1.0 / 3));
				cbrtXCbrtHiTab[i - MIN_DOUBLE_EXPONENT] = twoPowExpDiv3 * 0.5; // Half cbrt, to avoid overflows.
				cbrtSlopeHiTab[i - MIN_DOUBLE_EXPONENT] = (4.0 / 3) / (twoPowExpDiv3 * twoPowExpDiv3);
			}
			cbrtXCbrtLoTab[0] = 1.0;
			cbrtSlopeLoTab[0] = 1.0;
			final long CBRT_LO_MASK = (0x3FF0000000000000L | (0x000FFFFFFFFFFFFFL >> CBRT_LO_BITS));
			for (int i = 1; i < CBRT_LO_TAB_SIZE; i++) {
				long xBits = CBRT_LO_MASK | (((long) (i - 1)) << (52 - CBRT_LO_BITS));
				double cbrtX = StrictMath.cbrt(Double.longBitsToDouble(xBits));
				cbrtXCbrtLoTab[i] = cbrtX;
				cbrtSlopeLoTab[i] = 1 / (cbrtX * cbrtX);
			}
		}
	}

	// --------------------------------------------------------------------------
	// CONSTANTS FOR HYPOT
	// --------------------------------------------------------------------------

	/**
	 * For using sqrt, to avoid overflow/underflow, we want values magnitude in
	 * [1/sqrt(Double.MAX_VALUE/n),sqrt(Double.MAX_VALUE/n)], n being the number of
	 * arguments.
	 * 
	 * sqrt(Double.MAX_VALUE/2) = 9.480751908109176E153 and sqrt(Double.MAX_VALUE/3)
	 * = 7.741001517595157E153 so 2^511 = 6.7039039649712985E153 works for both.
	 */
	static final double HYPOT_MAX_MAG = NumbersUtils.twoPow(511);

	/**
	 * Large enough to get a value's magnitude back into [2^-511,2^511] from
	 * Double.MIN_VALUE or Double.MAX_VALUE, and small enough not to get it across
	 * that range (considering a 2*53 tolerance due to only checking magnitude of
	 * min/max value, and scaling all values together).
	 */
	static final double HYPOT_FACTOR = NumbersUtils.twoPow(750);

	// --------------------------------------------------------------------------
	// PUBLIC METHODS
	// --------------------------------------------------------------------------

	/**
	 * Ensures that all look-up tables are initialized - otherwise they are
	 * initialized lazily.
	 */
	public static void initTables() {
		/*
		 * Taking care not to call init methods here, which would recompute tables each
		 * time (even though the computations should be identical, since done with
		 * strictfp and StrictMath).
		 */
		int antiOptim = 0;
		antiOptim += MyTSinCos.sinTab.length;
		antiOptim += MyTTan.tanTab.length;
		antiOptim += MyTAsin.asinTab.length;
		antiOptim += MyTAsinPow.asinPowTab.length;
		antiOptim += MyTAtan.atanTab.length;
		antiOptim += MyTExp.expHiTab.length;
		antiOptim += MyTLog.logXTab.length;
		antiOptim += MyTTwoPow.twoPowTab.length;
		antiOptim += MyTSqrt.sqrtXSqrtHiTab.length;
		antiOptim += MyTCbrt.cbrtXCbrtHiTab.length;
		if (StrictMath.cos((double) antiOptim) == 0.0) {
			// Can't happen, cos is never +-0.0.
			throw new AssertionError();
		}
	}

	/*
	 * logarithms
	 */

	/**
	 * @param value An integer value in [1,Integer.MAX_VALUE].
	 * @return The integer part of the logarithm, in base 2, of the specified value,
	 *         i.e. a result in [0,30]
	 * @throws IllegalArgumentException if the specified value is <= 0.
	 */
	public static int log2(int value) {
		return NumbersUtils.log2(value);
	}

	/**
	 * @param value An integer value in [1,Long.MAX_VALUE].
	 * @return The integer part of the logarithm, in base 2, of the specified value,
	 *         i.e. a result in [0,62]
	 * @throws IllegalArgumentException if the specified value is <= 0.
	 */
	public static int log2(long value) {
		return NumbersUtils.log2(value);
	}

	/*
	 * powers
	 */

	/**
	 * Returns the exact result, provided it's in double range, i.e. if power is in
	 * [-1074,1023].
	 * 
	 * @param power An int power.
	 * @return 2^power as a double, or +-Infinity in case of overflow.
	 */
	public static double twoPow(int power) {
		/*
		 * OK to have this method factored here even though it returns a floating point
		 * value, because it only does integer operations and only takes integer
		 * arguments, so should behave the same even if inlined into FP-wide context.
		 */
		if (USE_TWO_POW_TAB) {
			if (power >= MIN_DOUBLE_EXPONENT) {
				if (power <= MAX_DOUBLE_EXPONENT) { // Normal or subnormal.
					return MyTTwoPow.twoPowTab[power - MIN_DOUBLE_EXPONENT];
				} else { // Overflow.
					return Double.POSITIVE_INFINITY;
				}
			} else { // Underflow.
				return 0.0;
			}
		} else {
			return NumbersUtils.twoPow(power);
		}
	}

	/**
	 * @param value An int value.
	 * @return value*value.
	 */
	public static int pow2(int value) {
		return value * value;
	}

	/**
	 * @param value A long value.
	 * @return value*value.
	 */
	public static long pow2(long value) {
		return value * value;
	}

	/**
	 * @param value An int value.
	 * @return value*value*value.
	 */
	public static int pow3(int value) {
		return value * value * value;
	}

	/**
	 * @param value A long value.
	 * @return value*value*value.
	 */
	public static long pow3(long value) {
		return value * value * value;
	}

	/*
	 * absolute values
	 */

	/**
	 * @param value An int value.
	 * @return The absolute value, except if value is Integer.MIN_VALUE, for which
	 *         it returns Integer.MIN_VALUE.
	 */
	public static int abs(int value) {
		if (FM_USE_JDK_MATH || SFM_USE_JDK_MATH) {
			return Math.abs(value);
		}
		return NumbersUtils.abs(value);
	}

	/**
	 * @param value A long value.
	 * @return The absolute value, except if value is Long.MIN_VALUE, for which it
	 *         returns Long.MIN_VALUE.
	 */
	public static long abs(long value) {
		if (FM_USE_JDK_MATH || SFM_USE_JDK_MATH) {
			return Math.abs(value);
		}
		return NumbersUtils.abs(value);
	}

	/*
	 * close values
	 */

	/**
	 * @param value A long value.
	 * @return The specified value as int.
	 * @throws ArithmeticException if the specified value is not in
	 *                             [Integer.MIN_VALUE,Integer.MAX_VALUE] range.
	 */
	public static int toIntExact(long value) {
		return NumbersUtils.asInt(value);
	}

	/**
	 * @param value A long value.
	 * @return The closest int value in [Integer.MIN_VALUE,Integer.MAX_VALUE] range.
	 */
	public static int toInt(long value) {
		return NumbersUtils.toInt(value);
	}

	/*
	 * ranges
	 */

	/**
	 * @param min   An int value.
	 * @param max   An int value.
	 * @param value An int value.
	 * @return minValue if value < minValue, maxValue if value > maxValue, value
	 *         otherwise.
	 */
	public static int toRange(int min, int max, int value) {
		return NumbersUtils.toRange(min, max, value);
	}

	/**
	 * @param min   A long value.
	 * @param max   A long value.
	 * @param value A long value.
	 * @return min if value < min, max if value > max, value otherwise.
	 */
	public static long toRange(long min, long max, long value) {
		return NumbersUtils.toRange(min, max, value);
	}

	/*
	 * unary operators (increment,decrement,negate)
	 */

	/**
	 * @param value An int value.
	 * @return The argument incremented by one.
	 * @throws ArithmeticException if the mathematical result is not in int range.
	 */
	public static int incrementExact(int value) {
		if (value == Integer.MAX_VALUE) {
			throw new ArithmeticException("integer overflow");
		}
		return value + 1;
	}

	/**
	 * @param value A long value.
	 * @return The argument incremented by one.
	 * @throws ArithmeticException if the mathematical result is not in long range.
	 */
	public static long incrementExact(long value) {
		if (value == Long.MAX_VALUE) {
			throw new ArithmeticException("long overflow");
		}
		return value + 1L;
	}

	/**
	 * @param value An int value.
	 * @return The argument incremented by one, or the argument if the mathematical
	 *         result is not in int range.
	 */
	public static int incrementBounded(int value) {
		if (value == Integer.MAX_VALUE) {
			return value;
		}
		return value + 1;
	}

	/**
	 * @param value A long value.
	 * @return The argument incremented by one, or the argument if the mathematical
	 *         result is not in long range.
	 */
	public static long incrementBounded(long value) {
		if (value == Long.MAX_VALUE) {
			return value;
		}
		return value + 1L;
	}

	/**
	 * @param value An int value.
	 * @return The argument decremented by one.
	 * @throws ArithmeticException if the mathematical result is not in int range.
	 */
	public static int decrementExact(int value) {
		if (value == Integer.MIN_VALUE) {
			throw new ArithmeticException("integer overflow");
		}
		return value - 1;
	}

	/**
	 * @param value A long value.
	 * @return The argument decremented by one.
	 * @throws ArithmeticException if the mathematical result is not in long range.
	 */
	public static long decrementExact(long value) {
		if (value == Long.MIN_VALUE) {
			throw new ArithmeticException("long overflow");
		}
		return value - 1L;
	}

	/**
	 * @param value An int value.
	 * @return The argument decremented by one, or the argument if the mathematical
	 *         result is not in int range.
	 */
	public static int decrementBounded(int value) {
		if (value == Integer.MIN_VALUE) {
			return value;
		}
		return value - 1;
	}

	/**
	 * @param value A long value.
	 * @return The argument decremented by one, or the argument if the mathematical
	 *         result is not in long range.
	 */
	public static long decrementBounded(long value) {
		if (value == Long.MIN_VALUE) {
			return value;
		}
		return value - 1L;
	}

	/**
	 * @param value An int value.
	 * @return The argument negated.
	 * @throws ArithmeticException if the mathematical result is not in int range.
	 */
	public static int negateExact(int value) {
		if (value == Integer.MIN_VALUE) {
			throw new ArithmeticException("integer overflow");
		}
		return -value;
	}

	/**
	 * @param value A long value.
	 * @return The argument negated.
	 * @throws ArithmeticException if the mathematical result is not in long range.
	 */
	public static long negateExact(long value) {
		if (value == Long.MIN_VALUE) {
			throw new ArithmeticException("long overflow");
		}
		return -value;
	}

	/**
	 * @param value An int value.
	 * @return The argument negated, or Integer.MAX_VALUE if the argument is
	 *         Integer.MIN_VALUE.
	 */
	public static int negateBounded(int value) {
		if (value == Integer.MIN_VALUE) {
			return Integer.MAX_VALUE;
		}
		return -value;
	}

	/**
	 * @param value A long value.
	 * @return The argument negated, or Long.MAX_VALUE if the argument is
	 *         Long.MIN_VALUE.
	 */
	public static long negateBounded(long value) {
		if (value == Long.MIN_VALUE) {
			return Long.MAX_VALUE;
		}
		return -value;
	}

	/*
	 * binary operators (+,-,*)
	 */

	/**
	 * @param a An int value.
	 * @param b An int value.
	 * @return The mathematical result of a+b.
	 * @throws ArithmeticException if the mathematical result of a+b is not in
	 *                             [Integer.MIN_VALUE,Integer.MAX_VALUE] range.
	 */
	public static int addExact(int a, int b) {
		return NumbersUtils.plusExact(a, b);
	}

	/**
	 * @param a A long value.
	 * @param b A long value.
	 * @return The mathematical result of a+b.
	 * @throws ArithmeticException if the mathematical result of a+b is not in
	 *                             [Long.MIN_VALUE,Long.MAX_VALUE] range.
	 */
	public static long addExact(long a, long b) {
		return NumbersUtils.plusExact(a, b);
	}

	/**
	 * @param a An int value.
	 * @param b An int value.
	 * @return The int value of [Integer.MIN_VALUE,Integer.MAX_VALUE] range which is
	 *         the closest to mathematical result of a+b.
	 */
	public static int addBounded(int a, int b) {
		return NumbersUtils.plusBounded(a, b);
	}

	/**
	 * @param a A long value.
	 * @param b A long value.
	 * @return The long value of [Long.MIN_VALUE,Long.MAX_VALUE] range which is the
	 *         closest to mathematical result of a+b.
	 */
	public static long addBounded(long a, long b) {
		return NumbersUtils.plusBounded(a, b);
	}

	/**
	 * @param a An int value.
	 * @param b An int value.
	 * @return The mathematical result of a-b.
	 * @throws ArithmeticException if the mathematical result of a-b is not in
	 *                             [Integer.MIN_VALUE,Integer.MAX_VALUE] range.
	 */
	public static int subtractExact(int a, int b) {
		return NumbersUtils.minusExact(a, b);
	}

	/**
	 * @param a A long value.
	 * @param b A long value.
	 * @return The mathematical result of a-b.
	 * @throws ArithmeticException if the mathematical result of a-b is not in
	 *                             [Long.MIN_VALUE,Long.MAX_VALUE] range.
	 */
	public static long subtractExact(long a, long b) {
		return NumbersUtils.minusExact(a, b);
	}

	/**
	 * @param a An int value.
	 * @param b An int value.
	 * @return The int value of [Integer.MIN_VALUE,Integer.MAX_VALUE] range which is
	 *         the closest to mathematical result of a-b.
	 */
	public static int subtractBounded(int a, int b) {
		return NumbersUtils.minusBounded(a, b);
	}

	/**
	 * @param a A long value.
	 * @param b A long value.
	 * @return The long value of [Long.MIN_VALUE,Long.MAX_VALUE] range which is the
	 *         closest to mathematical result of a-b.
	 */
	public static long subtractBounded(long a, long b) {
		return NumbersUtils.minusBounded(a, b);
	}

	/**
	 * @param a An int value.
	 * @param b An int value.
	 * @return The mathematical result of a*b.
	 * @throws ArithmeticException if the mathematical result of a*b is not in
	 *                             [Integer.MIN_VALUE,Integer.MAX_VALUE] range.
	 */
	public static int multiplyExact(int a, int b) {
		return NumbersUtils.timesExact(a, b);
	}

	/**
	 * @param a A long value.
	 * @param b An int value.
	 * @return The mathematical result of a*b.
	 * @throws ArithmeticException if the mathematical result of a*b is not in
	 *                             [Long.MIN_VALUE,Long.MAX_VALUE] range.
	 */
	public static long multiplyExact(long a, int b) {
		return NumbersUtils.timesExact(a, (long) b);
	}

	/**
	 * @param a A long value.
	 * @param b A long value.
	 * @return The mathematical result of a*b.
	 * @throws ArithmeticException if the mathematical result of a*b is not in
	 *                             [Long.MIN_VALUE,Long.MAX_VALUE] range.
	 */
	public static long multiplyExact(long a, long b) {
		return NumbersUtils.timesExact(a, b);
	}

	/**
	 * @param a An int value.
	 * @param b An int value.
	 * @return The int value of [Integer.MIN_VALUE,Integer.MAX_VALUE] range which is
	 *         the closest to mathematical result of a*b.
	 */
	public static int multiplyBounded(int a, int b) {
		return NumbersUtils.timesBounded(a, b);
	}

	/**
	 * @param a A long value.
	 * @param b An int value.
	 * @return The long value of [Long.MIN_VALUE,Long.MAX_VALUE] range which is the
	 *         closest to mathematical result of a*b.
	 */
	public static long multiplyBounded(long a, int b) {
		return NumbersUtils.timesBounded(a, (long) b);
	}

	/**
	 * @param a A long value.
	 * @param b A long value.
	 * @return The long value of [Long.MIN_VALUE,Long.MAX_VALUE] range which is the
	 *         closest to mathematical result of a*b.
	 */
	public static long multiplyBounded(long a, long b) {
		return NumbersUtils.timesBounded(a, b);
	}

	/**
	 * @param x An int value.
	 * @param y An int value.
	 * @return The mathematical product as a long.
	 */
	public static long multiplyFull(int x, int y) {
		return ((long) x) * ((long) y);
	}

	/**
	 * @param x A long value.
	 * @param y A long value.
	 * @return The most significant 64 bits of the 128-bit product of two 64-bit
	 *         factors.
	 */
	public static long multiplyHigh(long x, long y) {
		if ((x | y) < 0) {
			// Use technique from section 8-2 of Henry S. Warren, Jr.,
			// Hacker's Delight (2nd ed.) (Addison Wesley, 2013), 173-174.
			long x1 = (x >> 32);
			long y1 = (y >> 32);
			long x2 = (x & 0xFFFFFFFFL);
			long y2 = (y & 0xFFFFFFFFL);
			long z2 = x2 * y2;
			long t = x1 * y2 + (z2 >>> 32);
			long z1 = (t & 0xFFFFFFFFL) + x2 * y1;
			long z0 = (t >> 32);
			return x1 * y1 + z0 + (z1 >> 32);
		} else {
			// Use Karatsuba technique with two base 2^32 digits.
			long x1 = (x >>> 32);
			long y1 = (y >>> 32);
			long x2 = (x & 0xFFFFFFFFL);
			long y2 = (y & 0xFFFFFFFFL);
			long A = x1 * y1;
			long B = x2 * y2;
			long C = (x1 + x2) * (y1 + y2);
			long K = C - A - B;
			return (((B >>> 32) + K) >>> 32) + A;
		}
	}

	/*
	 * binary operators (/,%)
	 */

	/**
	 * Returns the largest int <= dividend/divisor.
	 * 
	 * Unlike "/" operator, which rounds towards 0, this division rounds towards
	 * -Infinity (which give different result when the exact result is negative).
	 * 
	 * @param x The dividend.
	 * @param y The divisor.
	 * @return The largest int <= dividend/divisor, unless dividend is
	 *         Integer.MIN_VALUE and divisor is -1, in which case Integer.MIN_VALUE
	 *         is returned.
	 * @throws ArithmeticException if the divisor is zero.
	 */
	public static int floorDiv(int x, int y) {
		int r = x / y;
		// If the signs are different and modulo not zero, rounding down.
		if (((x ^ y) < 0) && ((r * y) != x)) {
			r--;
		}
		return r;
	}

	/**
	 * Returns the largest long <= dividend/divisor.
	 * 
	 * Unlike "/" operator, which rounds towards 0, this division rounds towards
	 * -Infinity (which give different result when the exact result is negative).
	 * 
	 * @param x The dividend.
	 * @param y The divisor.
	 * @return The largest long <= dividend/divisor, unless dividend is
	 *         Long.MIN_VALUE and divisor is -1, in which case Long.MIN_VALUE is
	 *         returned.
	 * @throws ArithmeticException if the divisor is zero.
	 */
	public static long floorDiv(long x, int y) {
		return floorDiv(x, (long) y);
	}

	/**
	 * Returns the largest long <= dividend/divisor.
	 * 
	 * Unlike "/" operator, which rounds towards 0, this division rounds towards
	 * -Infinity (which give different result when the exact result is negative).
	 * 
	 * @param x The dividend.
	 * @param y The divisor.
	 * @return The largest long <= dividend/divisor, unless dividend is
	 *         Long.MIN_VALUE and divisor is -1, in which case Long.MIN_VALUE is
	 *         returned.
	 * @throws ArithmeticException if the divisor is zero.
	 */
	public static long floorDiv(long x, long y) {
		long r = x / y;
		// If the signs are different and modulo not zero, rounding down.
		if (((x ^ y) < 0) && ((r * y) != x)) {
			r--;
		}
		return r;
	}

	/**
	 * Returns the floor modulus, which is "x - floorDiv(x,y) * y", has the same
	 * sign as y, and is in ]-abs(y),abs(y)[.
	 *
	 * The relationship between floorMod and floorDiv is the same than between "%"
	 * and "/".
	 * 
	 * @param x The dividend.
	 * @param y The divisor.
	 * @return The floor modulus, i.e. "x - (floorDiv(x, y) * y)".
	 * @throws ArithmeticException if the divisor is zero.
	 */
	public static int floorMod(int x, int y) {
		return x - floorDiv(x, y) * y;
	}

	/**
	 * Returns the floor modulus, which is "x - floorDiv(x,y) * y", has the same
	 * sign as y, and is in ]-abs(y),abs(y)[.
	 *
	 * The relationship between floorMod and floorDiv is the same than between "%"
	 * and "/".
	 * 
	 * @param x The dividend.
	 * @param y The divisor.
	 * @return The floor modulus, i.e. "x - (floorDiv(x, y) * y)".
	 * @throws ArithmeticException if the divisor is zero.
	 */
	public static int floorMod(long x, int y) {
		// No overflow so can cast.
		return (int) (x - floorDiv(x, y) * y);
	}

	/**
	 * Returns the floor modulus, which is "x - floorDiv(x,y) * y", has the same
	 * sign as y, and is in ]-abs(y),abs(y)[.
	 *
	 * The relationship between floorMod and floorDiv is the same than between "%"
	 * and "/".
	 * 
	 * @param x The dividend.
	 * @param y The divisor.
	 * @return The floor modulus, i.e. "x - (floorDiv(x, y) * y)".
	 * @throws ArithmeticException if the divisor is zero.
	 */
	public static long floorMod(long x, long y) {
		return x - floorDiv(x, y) * y;
	}

	/*
	 * Non-redefined Math public values and treatments.
	 */

	public static int min(int a, int b) {
		return Math.min(a, b);
	}

	public static long min(long a, long b) {
		return Math.min(a, b);
	}

	public static int max(int a, int b) {
		return Math.max(a, b);
	}

	public static long max(long a, long b) {
		return Math.max(a, b);
	}

	// --------------------------------------------------------------------------
	// PACKAGE-PRIVATE METHODS
	// --------------------------------------------------------------------------

	/**
	 * @param power Must be in normal values range.
	 */
	static double twoPowNormal(int power) {
		if (USE_TWO_POW_TAB) {
			return MyTTwoPow.twoPowTab[power - MIN_DOUBLE_EXPONENT];
		} else {
			return Double.longBitsToDouble(((long) (power + MAX_DOUBLE_EXPONENT)) << 52);
		}
	}

	/**
	 * @param power Must be in normal or subnormal values range.
	 */
	static double twoPowNormalOrSubnormal(int power) {
		if (USE_TWO_POW_TAB) {
			return MyTTwoPow.twoPowTab[power - MIN_DOUBLE_EXPONENT];
		} else {
			if (power <= -MAX_DOUBLE_EXPONENT) { // Not normal.
				return Double.longBitsToDouble(0x0008000000000000L >> (-(power + MAX_DOUBLE_EXPONENT)));
			} else { // Normal.
				return Double.longBitsToDouble(((long) (power + MAX_DOUBLE_EXPONENT)) << 52);
			}
		}
	}

	static double atan2_pinf_yyy(double y) {
		if (y == Double.POSITIVE_INFINITY) {
			return Math.PI / 4;
		} else if (y == Double.NEGATIVE_INFINITY) {
			return -Math.PI / 4;
		} else if (y > 0.0) {
			return 0.0;
		} else if (y < 0.0) {
			return -0.0;
		} else {
			return Double.NaN;
		}
	}

	static double atan2_ninf_yyy(double y) {
		if (y == Double.POSITIVE_INFINITY) {
			return 3 * Math.PI / 4;
		} else if (y == Double.NEGATIVE_INFINITY) {
			return -3 * Math.PI / 4;
		} else if (y > 0.0) {
			return Math.PI;
		} else if (y < 0.0) {
			return -Math.PI;
		} else {
			return Double.NaN;
		}
	}

	static double atan2_yyy_zeroOrNaN(double y, double x) {
		if (x == 0.0) {
			if (y == 0.0) {
				if (signFromBit_antiCyclic(x) < 0) {
					// x is -0.0
					return signFromBit_antiCyclic(y) * Math.PI;
				} else {
					// +-0.0
					return y;
				}
			}
			if (y > 0.0) {
				return Math.PI / 2;
			} else if (y < 0.0) {
				return -Math.PI / 2;
			} else {
				return Double.NaN;
			}
		} else {
			return Double.NaN;
		}
	}

	/**
	 * At least one of the arguments must be NaN.
	 */
	static double hypot_NaN(double xAbs, double yAbs) {
		if ((xAbs == Double.POSITIVE_INFINITY) || (yAbs == Double.POSITIVE_INFINITY)) {
			return Double.POSITIVE_INFINITY;
		} else {
			return Double.NaN;
		}
	}

	/**
	 * At least one of the arguments must be NaN.
	 */
	static double hypot_NaN(double xAbs, double yAbs, double zAbs) {
		if ((xAbs == Double.POSITIVE_INFINITY) || (yAbs == Double.POSITIVE_INFINITY) || (zAbs == Double.POSITIVE_INFINITY)) {
			return Double.POSITIVE_INFINITY;
		} else {
			return Double.NaN;
		}
	}

	/*
	 * 
	 */

	/**
	 * @param remainder Must have 1 for 2nd and 3rd exponent bits, which is the case
	 *                  for heavyRemPiO2 remainders (their absolute values are >=
	 *                  Double.longBitsToDouble(0x3000000000000000L) =
	 *                  1.727233711018889E-77, and even if they were not, turning
	 *                  these bits from 0 to 1 on decoding would not change the
	 *                  absolute error much), and also works for +-Infinity or NaN
	 *                  encoding.
	 * @param quadrant  Must be in [0,3].
	 * @return Bits holding remainder, and quadrant instead of reamainder's 2nd and
	 *         3rd exponent bits.
	 */
	static long encodeRemainderAndQuadrant(double remainder, int quadrant) {
		final long bits = Double.doubleToRawLongBits(remainder);
		return (bits & QUADRANT_BITS_0_MASK) | (((long) quadrant) << 60);
	}

	static double decodeRemainder(long bits) {
		return Double.longBitsToDouble((bits & QUADRANT_BITS_0_MASK) | QUADRANT_PLACE_BITS);
	}

	static int decodeQuadrant(long bits) {
		return ((int) (bits >> 60)) & 3;
	}

	/*
	 * JDK-based remainders. Since a strict one for (% (PI/2)) is needed for
	 * heavyRemainderPiO2, we need it in this class. Then, for homogeneity, we put
	 * them all in this class. Then, to avoid code duplication for these slow-anyway
	 * methods, we just stick with strict versions, for both FastMath and
	 * StrictFastMath.
	 */

	/**
	 * @param angle Angle, in radians.
	 * @return Remainder of (angle % (2*PI)), in [-PI,PI].
	 */
	static strictfp double jdkRemainderTwoPi(double angle) {
		final double sin = StrictMath.sin(angle);
		final double cos = StrictMath.cos(angle);
		return StrictMath.atan2(sin, cos);
	}

	/**
	 * @param angle Angle, in radians.
	 * @return Remainder of (angle % PI), in [-PI/2,PI/2].
	 */
	static strictfp double jdkRemainderPi(double angle) {
		final double sin = StrictMath.sin(angle);
		final double cos = StrictMath.cos(angle);
		/*
		 * Making sure atan2's result ends up in [-PI/2,PI/2], i.e. has maximum
		 * accuracy.
		 */
		return StrictMath.atan2(sin, Math.abs(cos));
	}

	/**
	 * @param angle Angle, in radians.
	 * @return Bits of double corresponding to remainder of (angle % (PI/2)), in
	 *         [-PI/4,PI/4], with quadrant encoded in exponent bits.
	 */
	static strictfp long jdkRemainderPiO2(double angle, boolean negateRem) {
		final double sin = StrictMath.sin(angle);
		final double cos = StrictMath.cos(angle);

		/*
		 * Computing quadrant first, and then computing atan2, to make sure its result
		 * ends up in [-PI/4,PI/4], i.e. has maximum accuracy.
		 */

		final int q;
		final double sinForAtan2;
		final double cosForAtan2;
		if (cos >= (SQRT_2 / 2)) {
			// [-PI/4,PI/4]
			q = 0;
			sinForAtan2 = sin;
			cosForAtan2 = cos;
		} else if (cos <= -(SQRT_2 / 2)) {
			// [3*PI/4,5*PI/4]
			q = 2;
			sinForAtan2 = -sin;
			cosForAtan2 = -cos;
		} else if (sin > 0.0) {
			// [PI/4,3*PI/4]
			q = 1;
			sinForAtan2 = -cos;
			cosForAtan2 = sin;
		} else {
			// [5*PI/4,7*PI/4]
			q = 3;
			sinForAtan2 = cos;
			cosForAtan2 = -sin;
		}

		double fw = StrictMath.atan2(sinForAtan2, cosForAtan2);

		return encodeRemainderAndQuadrant(negateRem ? -fw : fw, q);
	}

	/*
	 * Our remainders implementations.
	 */

	/**
	 * @param angle Angle, in radians. Must not be NaN nor +-Infinity.
	 * @return Remainder of (angle % (2*PI)), in [-PI,PI].
	 */
	static strictfp double heavyRemainderTwoPi(double angle) {
		final long remAndQuad = heavyRemainderPiO2(angle, false);
		final double rem = decodeRemainder(remAndQuad);
		final int q = decodeQuadrant(remAndQuad);
		if (q == 0) {
			return rem;
		} else if (q == 1) {
			return (rem + PIO2_LO) + PIO2_HI;
		} else if (q == 2) {
			if (rem < 0.0) {
				return (rem + PI_LO) + PI_HI;
			} else {
				return (rem - PI_LO) - PI_HI;
			}
		} else {
			return (rem - PIO2_LO) - PIO2_HI;
		}
	}

	/**
	 * @param angle Angle, in radians. Must not be NaN nor +-Infinity.
	 * @return Remainder of (angle % PI), in [-PI/2,PI/2].
	 */
	static strictfp double heavyRemainderPi(double angle) {
		final long remAndQuad = heavyRemainderPiO2(angle, false);
		final double rem = decodeRemainder(remAndQuad);
		final int q = decodeQuadrant(remAndQuad);
		if ((q & 1) != 0) {
			// q is 1 or 3
			if (rem < 0.0) {
				return (rem + PIO2_LO) + PIO2_HI;
			} else {
				return (rem - PIO2_LO) - PIO2_HI;
			}
		}
		return rem;
	}

	/**
	 * Remainder using an accurate definition of PI. Derived from a fdlibm treatment
	 * called __kernel_rem_pio2.
	 * 
	 * Not defining a non-strictfp version for FastMath, to avoid duplicating its
	 * long and messy code, and because it's slow anyway, and should be rarely used
	 * when speed matters.
	 * 
	 * @param angle     Angle, in radians. Must not be NaN nor +-Infinity.
	 * @param negateRem True if remainder must be negated before encoded into
	 *                  returned long.
	 * @return Bits of double corresponding to remainder of (angle % (PI/2)), in
	 *         [-PI/4,PI/4], with quadrant encoded in exponent bits.
	 */
	static strictfp long heavyRemainderPiO2(double angle, boolean negateRem) {

		/*
		 * fdlibm treatments unrolled, to avoid garbage and be OOME-free, corresponding
		 * to: 1) initial jk = 4 (precision = 3 = 64 bits (extended)), which is more
		 * accurate than using precision = 2 (53 bits, double), even though we work with
		 * doubles and use strictfp! 2) max lengths of 8 for f[], 6 for q[], fq[] and
		 * iq[]. 3) at most one recomputation (one goto). These limitations were
		 * experimentally found to be sufficient for billions of random doubles of
		 * random magnitudes. For the rare cases that our unrolled treatments can't
		 * handle, we fall back to a JDK-based implementation.
		 */

		int n, i, j, ih;
		double fw;

		/*
		 * Turning angle into 24-bits integer chunks. Done outside __kernel_rem_pio2,
		 * but we factor it inside our method.
		 */

		// Reworking exponent to have a value < 2^24.
		final long lx = Double.doubleToRawLongBits(angle);
		final long exp = ((lx >> 52) & 0x7FF) - (1023 + 23);
		double z = Double.longBitsToDouble(lx - (exp << 52));

		double x0 = (double) (int) z;
		z = (z - x0) * TWO_POW_24;
		double x1 = (double) (int) z;
		z = (z - x1) * TWO_POW_24;
		double x2 = (double) (int) z;

		final int e0 = (int) exp;
		// in [1,3]
		final int nx = (x2 == 0.0) ? ((x1 == 0.0) ? 1 : 2) : 3;

		/*
		 * 
		 */

		double f0, f1, f2, f3, f4, f5, f6, f7;
		double q0, q1, q2, q3, q4, q5;
		int iq0, iq1, iq2, iq3, iq4, iq5;

		int jk = 4;

		int jx = nx - 1;
		int jv = Math.max(0, (e0 - 3) / 24);
		// In fdlibm, this is q0, but we prefer to use q0 for q[0].
		int qZero = e0 - 24 * (jv + 1);

		j = jv - jx;
		if (jx == 0) {
			f6 = 0.0;
			f5 = 0.0;
			f4 = (j >= -4) ? TWO_OVER_PI_TAB[j + 4] : 0.0;
			f3 = (j >= -3) ? TWO_OVER_PI_TAB[j + 3] : 0.0;
			f2 = (j >= -2) ? TWO_OVER_PI_TAB[j + 2] : 0.0;
			f1 = (j >= -1) ? TWO_OVER_PI_TAB[j + 1] : 0.0;
			f0 = (j >= 0) ? TWO_OVER_PI_TAB[j] : 0.0;

			q0 = x0 * f0;
			q1 = x0 * f1;
			q2 = x0 * f2;
			q3 = x0 * f3;
			q4 = x0 * f4;
		} else if (jx == 1) {
			f6 = 0.0;
			f5 = (j >= -5) ? TWO_OVER_PI_TAB[j + 5] : 0.0;
			f4 = (j >= -4) ? TWO_OVER_PI_TAB[j + 4] : 0.0;
			f3 = (j >= -3) ? TWO_OVER_PI_TAB[j + 3] : 0.0;
			f2 = (j >= -2) ? TWO_OVER_PI_TAB[j + 2] : 0.0;
			f1 = (j >= -1) ? TWO_OVER_PI_TAB[j + 1] : 0.0;
			f0 = (j >= 0) ? TWO_OVER_PI_TAB[j] : 0.0;

			q0 = x0 * f1 + x1 * f0;
			q1 = x0 * f2 + x1 * f1;
			q2 = x0 * f3 + x1 * f2;
			q3 = x0 * f4 + x1 * f3;
			q4 = x0 * f5 + x1 * f4;
		} else { // jx == 2
			f6 = (j >= -6) ? TWO_OVER_PI_TAB[j + 6] : 0.0;
			f5 = (j >= -5) ? TWO_OVER_PI_TAB[j + 5] : 0.0;
			f4 = (j >= -4) ? TWO_OVER_PI_TAB[j + 4] : 0.0;
			f3 = (j >= -3) ? TWO_OVER_PI_TAB[j + 3] : 0.0;
			f2 = (j >= -2) ? TWO_OVER_PI_TAB[j + 2] : 0.0;
			f1 = (j >= -1) ? TWO_OVER_PI_TAB[j + 1] : 0.0;
			f0 = (j >= 0) ? TWO_OVER_PI_TAB[j] : 0.0;

			q0 = x0 * f2 + x1 * f1 + x2 * f0;
			q1 = x0 * f3 + x1 * f2 + x2 * f1;
			q2 = x0 * f4 + x1 * f3 + x2 * f2;
			q3 = x0 * f5 + x1 * f4 + x2 * f3;
			q4 = x0 * f6 + x1 * f5 + x2 * f4;
		}

		double twoPowQZero = twoPowNormal(qZero);

		int jz = jk;

		/*
		 * Unrolling of first round.
		 */

		z = q4;
		fw = (double) (int) (TWO_POW_N24 * z);
		iq0 = (int) (z - TWO_POW_24 * fw);
		z = q3 + fw;
		fw = (double) (int) (TWO_POW_N24 * z);
		iq1 = (int) (z - TWO_POW_24 * fw);
		z = q2 + fw;
		fw = (double) (int) (TWO_POW_N24 * z);
		iq2 = (int) (z - TWO_POW_24 * fw);
		z = q1 + fw;
		fw = (double) (int) (TWO_POW_N24 * z);
		iq3 = (int) (z - TWO_POW_24 * fw);
		z = q0 + fw;
		iq4 = 0;
		iq5 = 0;

		z = (z * twoPowQZero) % 8.0;
		n = (int) z;
		z -= (double) n;

		ih = 0;
		if (qZero > 0) {
			// Parentheses against code formatter bug.
			i = (iq3 >> (24 - qZero));
			n += i;
			iq3 -= i << (24 - qZero);
			ih = iq3 >> (23 - qZero);
		} else if (qZero == 0) {
			ih = iq3 >> 23;
		} else if (z >= 0.5) {
			ih = 2;
		}

		if (ih > 0) {
			n += 1;
			// carry = 1 is common case,
			// so using it as initial value.
			int carry = 1;
			if (iq0 != 0) {
				iq0 = 0x1000000 - iq0;
				iq1 = 0xFFFFFF - iq1;
				iq2 = 0xFFFFFF - iq2;
				iq3 = 0xFFFFFF - iq3;
			} else if (iq1 != 0) {
				iq1 = 0x1000000 - iq1;
				iq2 = 0xFFFFFF - iq2;
				iq3 = 0xFFFFFF - iq3;
			} else if (iq2 != 0) {
				iq2 = 0x1000000 - iq2;
				iq3 = 0xFFFFFF - iq3;
			} else if (iq3 != 0) {
				iq3 = 0x1000000 - iq3;
			} else {
				carry = 0;
			}
			if (qZero > 0) {
				if (qZero == 1) {
					iq3 &= 0x7FFFFF;
				} else if (qZero == 2) {
					iq3 &= 0x3FFFFF;
				}
			}
			if (ih == 2) {
				z = 1.0 - z;
				if (carry != 0) {
					z -= twoPowQZero;
				}
			}
		}

		if (z == 0.0) {
			if (iq3 == 0) {
				// With random values of random magnitude,
				// probability for this to happen seems lower than 1e-6.
				// jz would be more than just incremented by one,
				// which our unrolling doesn't support.
				return jdkRemainderPiO2(angle, negateRem);
			}
			if (jx == 0) {
				f5 = TWO_OVER_PI_TAB[jv + 5];
				q5 = x0 * f5;
			} else if (jx == 1) {
				f6 = TWO_OVER_PI_TAB[jv + 5];
				q5 = x0 * f6 + x1 * f5;
			} else { // jx == 2
				f7 = TWO_OVER_PI_TAB[jv + 5];
				q5 = x0 * f7 + x1 * f6 + x2 * f5;
			}

			jz++;

			/*
			 * Unrolling of second round.
			 */

			z = q5;
			fw = (double) (int) (TWO_POW_N24 * z);
			iq0 = (int) (z - TWO_POW_24 * fw);
			z = q4 + fw;
			fw = (double) (int) (TWO_POW_N24 * z);
			iq1 = (int) (z - TWO_POW_24 * fw);
			z = q3 + fw;
			fw = (double) (int) (TWO_POW_N24 * z);
			iq2 = (int) (z - TWO_POW_24 * fw);
			z = q2 + fw;
			fw = (double) (int) (TWO_POW_N24 * z);
			iq3 = (int) (z - TWO_POW_24 * fw);
			z = q1 + fw;
			fw = (double) (int) (TWO_POW_N24 * z);
			iq4 = (int) (z - TWO_POW_24 * fw);
			z = q0 + fw;
			iq5 = 0;

			z = (z * twoPowQZero) % 8.0;
			n = (int) z;
			z -= (double) n;

			ih = 0;
			if (qZero > 0) {
				// Parentheses against code formatter bug.
				i = (iq4 >> (24 - qZero));
				n += i;
				iq4 -= i << (24 - qZero);
				ih = iq4 >> (23 - qZero);
			} else if (qZero == 0) {
				ih = iq4 >> 23;
			} else if (z >= 0.5) {
				ih = 2;
			}

			if (ih > 0) {
				n += 1;
				// carry = 1 is common case,
				// so using it as initial value.
				int carry = 1;
				if (iq0 != 0) {
					iq0 = 0x1000000 - iq0;
					iq1 = 0xFFFFFF - iq1;
					iq2 = 0xFFFFFF - iq2;
					iq3 = 0xFFFFFF - iq3;
					iq4 = 0xFFFFFF - iq4;
				} else if (iq1 != 0) {
					iq1 = 0x1000000 - iq1;
					iq2 = 0xFFFFFF - iq2;
					iq3 = 0xFFFFFF - iq3;
					iq4 = 0xFFFFFF - iq4;
				} else if (iq2 != 0) {
					iq2 = 0x1000000 - iq2;
					iq3 = 0xFFFFFF - iq3;
					iq4 = 0xFFFFFF - iq4;
				} else if (iq3 != 0) {
					iq3 = 0x1000000 - iq3;
					iq4 = 0xFFFFFF - iq4;
				} else if (iq4 != 0) {
					iq4 = 0x1000000 - iq4;
				} else {
					carry = 0;
				}
				if (qZero > 0) {
					if (qZero == 1) {
						iq4 &= 0x7FFFFF;
					} else if (qZero == 2) {
						iq4 &= 0x3FFFFF;
					}
				}
				if (ih == 2) {
					z = 1.0 - z;
					if (carry != 0) {
						z -= twoPowQZero;
					}
				}
			}

			if (z == 0.0) {
				if (iq4 == 0) {
					// Case not encountered in tests, but still handling it.
					// Would require a third loop unrolling.
					return jdkRemainderPiO2(angle, negateRem);
				} else {
					// z == 0.0, and iq4 != 0,
					// so we remove 24 from qZero only once,
					// but since we no longer use qZero,
					// we just bother to multiply its 2-power
					// by 2^-24.
					jz--;
					twoPowQZero *= TWO_POW_N24;
				}
			} else {
				// z != 0.0 at end of second round.
			}
		} else {
			// z != 0.0 at end of first round.
		}

		/*
		 * After loop.
		 */

		if (z != 0.0) {
			z /= twoPowQZero;
			if (z >= TWO_POW_24) {
				fw = (double) (int) (TWO_POW_N24 * z);
				if (jz == jk) {
					iq4 = (int) (z - TWO_POW_24 * fw);
					jz++; // jz to 5
					// Not using qZero anymore so not updating it.
					twoPowQZero *= TWO_POW_24;
					iq5 = (int) fw;
				} else { // jz == jk+1 == 5
					// Case not encountered in tests, but still handling it.
					// Would require use of iq6, with jz = 6.
					return jdkRemainderPiO2(angle, negateRem);
				}
			} else {
				if (jz == jk) {
					iq4 = (int) z;
				} else { // jz == jk+1 == 5
					// Case not encountered in tests, but still handling it.
					iq5 = (int) z;
				}
			}
		}

		fw = twoPowQZero;

		if (jz == 5) {
			q5 = fw * (double) iq5;
			fw *= TWO_POW_N24;
		} else {
			q5 = 0.0;
		}
		q4 = fw * (double) iq4;
		fw *= TWO_POW_N24;
		q3 = fw * (double) iq3;
		fw *= TWO_POW_N24;
		q2 = fw * (double) iq2;
		fw *= TWO_POW_N24;
		q1 = fw * (double) iq1;
		fw *= TWO_POW_N24;
		q0 = fw * (double) iq0;

		/*
		 * We just use HI part of the result.
		 */

		fw = PIO2_TAB0 * q5;
		fw += PIO2_TAB0 * q4 + PIO2_TAB1 * q5;
		fw += PIO2_TAB0 * q3 + PIO2_TAB1 * q4 + PIO2_TAB2 * q5;
		fw += PIO2_TAB0 * q2 + PIO2_TAB1 * q3 + PIO2_TAB2 * q4 + PIO2_TAB3 * q5;
		fw += PIO2_TAB0 * q1 + PIO2_TAB1 * q2 + PIO2_TAB2 * q3 + PIO2_TAB3 * q4 + PIO2_TAB4 * q5;
		fw += PIO2_TAB0 * q0 + PIO2_TAB1 * q1 + PIO2_TAB2 * q2 + PIO2_TAB3 * q3 + PIO2_TAB4 * q4 + PIO2_TAB5 * q5;

		if ((ih != 0) ^ negateRem) {
			fw = -fw;
		}

		return encodeRemainderAndQuadrant(fw, n & 3);
	}

	// --------------------------------------------------------------------------
	// PRIVATE METHODS
	// --------------------------------------------------------------------------

	/**
	 * Redefined here, to avoid cyclic dependency with (Strict)FastMath.
	 * 
	 * @param value A double value.
	 * @return -1 if sign bit is 1, 1 if sign bit is 0.
	 */
	private static long signFromBit_antiCyclic(double value) {
		// Returning a long, to avoid useless cast into int.
		return ((Double.doubleToRawLongBits(value) >> 62) | 1);
	}

	private static boolean getBooleanProperty(final String key, boolean defaultValue) {
		final String tmp = System.getProperty(key);
		if (tmp != null) {
			return Boolean.parseBoolean(tmp);
		} else {
			return defaultValue;
		}
	}

	/**
	 * Use look-up tables size power through this method, to make sure is it small
	 * in case java.lang.Math is directly used.
	 */
	private static int getTabSizePower(int tabSizePower) {
		return (FM_USE_JDK_MATH && SFM_USE_JDK_MATH) ? Math.min(2, tabSizePower) : tabSizePower;
	}

	// --------------------------------------------------------------------------
	// CONFIGURATION
	// --------------------------------------------------------------------------

	private static final boolean USE_JDK_MATH = FM_USE_JDK_MATH;

	private static final boolean USE_REDEFINED_LOG = FM_USE_REDEFINED_LOG;

	private static final boolean USE_REDEFINED_SQRT = FM_USE_REDEFINED_SQRT;

	private static final boolean USE_POWTABS_FOR_ASIN = FM_USE_POWTABS_FOR_ASIN;

	// --------------------------------------------------------------------------
	// PUBLIC METHODS
	// --------------------------------------------------------------------------

	/*
	 * trigonometry
	 */

	/**
	 * @param angle Angle in radians.
	 * @return Angle sine.
	 */
	public static double sin(double angle) {
		if (USE_JDK_MATH) {
			return Math.sin(angle);
		}
		boolean negateResult = false;
		if (angle < 0.0) {
			angle = -angle;
			negateResult = true;
		}
		if (angle > SIN_COS_MAX_VALUE_FOR_INT_MODULO) {
			if (false) {
				// Can give very bad relative error near PI (mod 2*PI).
				angle = remainderTwoPi(angle);
				if (angle < 0.0) {
					angle = -angle;
					negateResult = !negateResult;
				}
			} else {
				final long remAndQuad = remainderPiO2(angle);
				angle = decodeRemainder(remAndQuad);
				final double sin;
				final int q = decodeQuadrant(remAndQuad);
				if (q == 0) {
					sin = sin(angle);
				} else if (q == 1) {
					sin = cos(angle);
				} else if (q == 2) {
					sin = -sin(angle);
				} else {
					sin = -cos(angle);
				}
				return (negateResult ? -sin : sin);
			}
		}
		// index: possibly outside tables range.
		int index = (int) (angle * SIN_COS_INDEXER + 0.5);
		double delta = (angle - index * SIN_COS_DELTA_HI) - index * SIN_COS_DELTA_LO;
		// Making sure index is within tables range.
		// Last value of each table is the same than first,
		// so we ignore it (tabs size minus one) for modulo.
		index &= (SIN_COS_TABS_SIZE - 2); // index % (SIN_COS_TABS_SIZE-1)
		double indexSin = MyTSinCos.sinTab[index];
		double indexCos = MyTSinCos.cosTab[index];
		double result = indexSin + delta * (indexCos + delta * (-indexSin * ONE_DIV_F2 + delta * (-indexCos * ONE_DIV_F3 + delta * indexSin * ONE_DIV_F4)));
		return negateResult ? -result : result;
	}

	/**
	 * Quick sin, with accuracy of about 1.6e-3 (PI/<look-up tabs size>) for |angle|
	 * < 6588395.0 (Integer.MAX_VALUE * (2*PI/<look-up tabs size>) - 2) (- 2 due to
	 * removing PI/2 before using cosine tab), and no accuracy at all for larger
	 * values.
	 * 
	 * @param angle Angle in radians.
	 * @return Angle sine.
	 */
	public static double sinQuick(double angle) {
		if (USE_JDK_MATH) {
			return Math.sin(angle);
		}
		return MyTSinCos.cosTab[((int) (Math.abs(angle - Math.PI / 2) * SIN_COS_INDEXER + 0.5)) & (SIN_COS_TABS_SIZE - 2)];
	}

	/**
	 * @param angle Angle in radians.
	 * @return Angle cosine.
	 */
	public static double cos(double angle) {
		if (USE_JDK_MATH) {
			return Math.cos(angle);
		}
		angle = Math.abs(angle);
		if (angle > SIN_COS_MAX_VALUE_FOR_INT_MODULO) {
			if (false) {
				// Can give very bad relative error near PI (mod 2*PI).
				angle = remainderTwoPi(angle);
				if (angle < 0.0) {
					angle = -angle;
				}
			} else {
				final long remAndQuad = remainderPiO2(angle);
				angle = decodeRemainder(remAndQuad);
				final double cos;
				final int q = decodeQuadrant(remAndQuad);
				if (q == 0) {
					cos = cos(angle);
				} else if (q == 1) {
					cos = -sin(angle);
				} else if (q == 2) {
					cos = -cos(angle);
				} else {
					cos = sin(angle);
				}
				return cos;
			}
		}
		// index: possibly outside tables range.
		int index = (int) (angle * SIN_COS_INDEXER + 0.5);
		double delta = (angle - index * SIN_COS_DELTA_HI) - index * SIN_COS_DELTA_LO;
		// Making sure index is within tables range.
		// Last value of each table is the same than first,
		// so we ignore it (tabs size minus one) for modulo.
		index &= (SIN_COS_TABS_SIZE - 2); // index % (SIN_COS_TABS_SIZE-1)
		double indexCos = MyTSinCos.cosTab[index];
		double indexSin = MyTSinCos.sinTab[index];
		return indexCos + delta * (-indexSin + delta * (-indexCos * ONE_DIV_F2 + delta * (indexSin * ONE_DIV_F3 + delta * indexCos * ONE_DIV_F4)));
	}

	/**
	 * Quick cos, with accuracy of about 1.6e-3 (PI/<look-up tabs size>) for |angle|
	 * < 6588397.0 (Integer.MAX_VALUE * (2*PI/<look-up tabs size>)), and no accuracy
	 * at all for larger values.
	 * 
	 * @param angle Angle in radians.
	 * @return Angle cosine.
	 */
	public static double cosQuick(double angle) {
		if (USE_JDK_MATH) {
			return Math.cos(angle);
		}
		return MyTSinCos.cosTab[((int) (Math.abs(angle) * SIN_COS_INDEXER + 0.5)) & (SIN_COS_TABS_SIZE - 2)];
	}

	/**
	 * Computes sine and cosine together.
	 * 
	 * @param angle  Angle in radians.
	 * @param cosine (out) Angle cosine.
	 * @return Angle sine.
	 */
	public static double sinAndCos(double angle, DoubleWrapper cosine) {
		if (USE_JDK_MATH) {
			cosine.value = Math.cos(angle);
			return Math.sin(angle);
		}
		// Using the same algorithm than sin(double) method,
		// and computing also cosine at the end.
		boolean negateResult = false;
		if (angle < 0.0) {
			angle = -angle;
			negateResult = true;
		}
		if (angle > SIN_COS_MAX_VALUE_FOR_INT_MODULO) {
			if (false) {
				// Can give very bad relative error near PI (mod 2*PI).
				angle = remainderTwoPi(angle);
				if (angle < 0.0) {
					angle = -angle;
					negateResult = !negateResult;
				}
			} else {
				final long remAndQuad = remainderPiO2(angle);
				angle = decodeRemainder(remAndQuad);
				final double sin;
				final int q = decodeQuadrant(remAndQuad);
				if (q == 0) {
					sin = sin(angle);
					cosine.value = cos(angle);
				} else if (q == 1) {
					sin = cos(angle);
					cosine.value = -sin(angle);
				} else if (q == 2) {
					sin = -sin(angle);
					cosine.value = -cos(angle);
				} else {
					sin = -cos(angle);
					cosine.value = sin(angle);
				}
				return (negateResult ? -sin : sin);
			}
		}
		int index = (int) (angle * SIN_COS_INDEXER + 0.5);
		double delta = (angle - index * SIN_COS_DELTA_HI) - index * SIN_COS_DELTA_LO;
		index &= (SIN_COS_TABS_SIZE - 2); // index % (SIN_COS_TABS_SIZE-1)
		double indexSin = MyTSinCos.sinTab[index];
		double indexCos = MyTSinCos.cosTab[index];
		// Could factor some multiplications (delta * factorials), but then is less
		// accurate.
		cosine.value = indexCos + delta * (-indexSin + delta * (-indexCos * ONE_DIV_F2 + delta * (indexSin * ONE_DIV_F3 + delta * indexCos * ONE_DIV_F4)));
		double result = indexSin + delta * (indexCos + delta * (-indexSin * ONE_DIV_F2 + delta * (-indexCos * ONE_DIV_F3 + delta * indexSin * ONE_DIV_F4)));
		return negateResult ? -result : result;
	}

	/**
	 * Can have very bad relative error near +-PI/2, but of the same magnitude than
	 * the relative delta between StrictMath.tan(PI/2) and
	 * StrictMath.tan(nextDown(PI/2)).
	 * 
	 * @param angle Angle in radians.
	 * @return Angle tangent.
	 */
	public static double tan(double angle) {
		if (USE_JDK_MATH) {
			return Math.tan(angle);
		}
		boolean negateResult = false;
		if (angle < 0.0) {
			angle = -angle;
			negateResult = true;
		}
		if (angle > TAN_MAX_VALUE_FOR_INT_MODULO) {
			angle = remainderPi(angle);
			if (angle < 0.0) {
				angle = -angle;
				negateResult = !negateResult;
			}
		}
		// index: possibly outside tables range.
		int index = (int) (angle * TAN_INDEXER + 0.5);
		double delta = (angle - index * TAN_DELTA_HI) - index * TAN_DELTA_LO;
		// Making sure index is within tables range.
		// index modulo PI, i.e. 2*(virtual tab size minus one).
		index &= (2 * (TAN_VIRTUAL_TABS_SIZE - 1) - 1); // index % (2*(TAN_VIRTUAL_TABS_SIZE-1))
		// Here, index is in [0,2*(TAN_VIRTUAL_TABS_SIZE-1)-1], i.e. indicates an angle
		// in [0,PI[.
		if (index > (TAN_VIRTUAL_TABS_SIZE - 1)) {
			index = (2 * (TAN_VIRTUAL_TABS_SIZE - 1)) - index;
			delta = -delta;
			negateResult = !negateResult;
		}
		double result;
		if (index < TAN_TABS_SIZE) {
			result = MyTTan.tanTab[index] + delta * (MyTTan.tanDer1DivF1Tab[index]
					+ delta * (MyTTan.tanDer2DivF2Tab[index] + delta * (MyTTan.tanDer3DivF3Tab[index] + delta * MyTTan.tanDer4DivF4Tab[index])));
		} else { // angle in ]TAN_MAX_VALUE_FOR_TABS,TAN_MAX_VALUE_FOR_INT_MODULO], or angle is
					// NaN
			// Using tan(angle) == 1/tan(PI/2-angle) formula: changing angle (index and
			// delta), and inverting.
			index = (TAN_VIRTUAL_TABS_SIZE - 1) - index;
			result = 1 / (MyTTan.tanTab[index] - delta * (MyTTan.tanDer1DivF1Tab[index]
					- delta * (MyTTan.tanDer2DivF2Tab[index] - delta * (MyTTan.tanDer3DivF3Tab[index] - delta * MyTTan.tanDer4DivF4Tab[index]))));
		}
		return negateResult ? -result : result;
	}

	/**
	 * @param value Value in [-1,1].
	 * @return Value arcsine, in radians, in [-PI/2,PI/2].
	 */
	public static double asin(double value) {
		if (USE_JDK_MATH) {
			return Math.asin(value);
		}
		boolean negateResult = false;
		if (value < 0.0) {
			value = -value;
			negateResult = true;
		}
		if (value <= ASIN_MAX_VALUE_FOR_TABS) {
			int index = (int) (value * ASIN_INDEXER + 0.5);
			double delta = value - index * ASIN_DELTA;
			double result = MyTAsin.asinTab[index] + delta * (MyTAsin.asinDer1DivF1Tab[index]
					+ delta * (MyTAsin.asinDer2DivF2Tab[index] + delta * (MyTAsin.asinDer3DivF3Tab[index] + delta * MyTAsin.asinDer4DivF4Tab[index])));
			return negateResult ? -result : result;
		} else if (USE_POWTABS_FOR_ASIN && (value <= ASIN_MAX_VALUE_FOR_POWTABS)) {
			int index = (int) (powFast(value * ASIN_POWTABS_ONE_DIV_MAX_VALUE, ASIN_POWTABS_POWER) * ASIN_POWTABS_SIZE_MINUS_ONE + 0.5);
			double delta = value - MyTAsinPow.asinParamPowTab[index];
			double result = MyTAsinPow.asinPowTab[index] + delta * (MyTAsinPow.asinDer1DivF1PowTab[index] + delta * (MyTAsinPow.asinDer2DivF2PowTab[index]
					+ delta * (MyTAsinPow.asinDer3DivF3PowTab[index] + delta * MyTAsinPow.asinDer4DivF4PowTab[index])));
			return negateResult ? -result : result;
		} else { // value > ASIN_MAX_VALUE_FOR_TABS, or value is NaN
			// This part is derived from fdlibm.
			if (value < 1.0) {
				double t = (1.0 - value) * 0.5;
				double p = t * (ASIN_PS0 + t * (ASIN_PS1 + t * (ASIN_PS2 + t * (ASIN_PS3 + t * (ASIN_PS4 + t * ASIN_PS5)))));
				double q = 1.0 + t * (ASIN_QS1 + t * (ASIN_QS2 + t * (ASIN_QS3 + t * ASIN_QS4)));
				double s = sqrt(t);
				double z = s + s * (p / q);
				double result = ASIN_PIO2_HI - ((z + z) - ASIN_PIO2_LO);
				return negateResult ? -result : result;
			} else { // value >= 1.0, or value is NaN
				if (value == 1.0) {
					return negateResult ? -Math.PI / 2 : Math.PI / 2;
				} else {
					return Double.NaN;
				}
			}
		}
	}

	/**
	 * If value is not NaN and is outside [-1,1] range, closest value in this range
	 * is used.
	 * 
	 * @param value Value in [-1,1].
	 * @return Value arcsine, in radians, in [-PI/2,PI/2].
	 */
	public static double asinInRange(double value) {
		if (value <= -1.0) {
			return -Math.PI / 2;
		} else if (value >= 1.0) {
			return Math.PI / 2;
		} else {
			return asin(value);
		}
	}

	/**
	 * @param value Value in [-1,1].
	 * @return Value arccosine, in radians, in [0,PI].
	 */
	public static double acos(double value) {
		if (USE_JDK_MATH) {
			return Math.acos(value);
		}
		return Math.PI / 2 - asin(value);
	}

	/**
	 * If value is not NaN and is outside [-1,1] range, closest value in this range
	 * is used.
	 * 
	 * @param value Value in [-1,1].
	 * @return Value arccosine, in radians, in [0,PI].
	 */
	public static double acosInRange(double value) {
		if (value <= -1.0) {
			return Math.PI;
		} else if (value >= 1.0) {
			return 0.0;
		} else {
			return acos(value);
		}
	}

	/**
	 * @param value A double value.
	 * @return Value arctangent, in radians, in [-PI/2,PI/2].
	 */
	public static double atan(double value) {
		if (USE_JDK_MATH) {
			return Math.atan(value);
		}
		boolean negateResult = false;
		if (value < 0.0) {
			value = -value;
			negateResult = true;
		}
		if (value == 1.0) {
			// We want "exact" result for 1.0.
			return negateResult ? -Math.PI / 4 : Math.PI / 4;
		} else if (value <= ATAN_MAX_VALUE_FOR_TABS) {
			int index = (int) (value * ATAN_INDEXER + 0.5);
			double delta = value - index * ATAN_DELTA;
			double result = MyTAtan.atanTab[index] + delta * (MyTAtan.atanDer1DivF1Tab[index]
					+ delta * (MyTAtan.atanDer2DivF2Tab[index] + delta * (MyTAtan.atanDer3DivF3Tab[index] + delta * MyTAtan.atanDer4DivF4Tab[index])));
			return negateResult ? -result : result;
		} else { // value > ATAN_MAX_VALUE_FOR_TABS, or value is NaN
			// This part is derived from fdlibm.
			if (value < TWO_POW_66) {
				double x = -1 / value;
				double x2 = x * x;
				double x4 = x2 * x2;
				double s1 = x2 * (ATAN_AT0 + x4 * (ATAN_AT2 + x4 * (ATAN_AT4 + x4 * (ATAN_AT6 + x4 * (ATAN_AT8 + x4 * ATAN_AT10)))));
				double s2 = x4 * (ATAN_AT1 + x4 * (ATAN_AT3 + x4 * (ATAN_AT5 + x4 * (ATAN_AT7 + x4 * ATAN_AT9))));
				double result = ATAN_HI3 - ((x * (s1 + s2) - ATAN_LO3) - x);
				return negateResult ? -result : result;
			} else { // value >= 2^66, or value is NaN
				if (value != value) {
					return Double.NaN;
				} else {
					return negateResult ? -Math.PI / 2 : Math.PI / 2;
				}
			}
		}
	}

	/**
	 * For special values for which multiple conventions could be adopted, behaves
	 * like Math.atan2(double,double).
	 * 
	 * @param y Coordinate on y axis.
	 * @param x Coordinate on x axis.
	 * @return Angle from x axis positive side to (x,y) position, in radians, in
	 *         [-PI,PI]. Angle measure is positive when going from x axis to y axis
	 *         (positive sides).
	 */
	public static double atan2(double y, double x) {
		if (USE_JDK_MATH) {
			return Math.atan2(y, x);
		}
		/*
		 * Using sub-methods, to make method lighter for general case, and to avoid
		 * JIT-optimization crash on NaN.
		 */
		if (x > 0.0) {
			if (y == 0.0) {
				// +-0.0
				return y;
			}
			if (x == Double.POSITIVE_INFINITY) {
				return atan2_pinf_yyy(y);
			} else {
				return atan(y / x);
			}
		} else if (x < 0.0) {
			if (y == 0.0) {
				return signFromBit(y) * Math.PI;
			}
			if (x == Double.NEGATIVE_INFINITY) {
				return atan2_ninf_yyy(y);
			} else if (y > 0.0) {
				return Math.PI / 2 - atan(x / y);
			} else if (y < 0.0) {
				return -Math.PI / 2 - atan(x / y);
			} else {
				return Double.NaN;
			}
		} else {
			return atan2_yyy_zeroOrNaN(y, x);
		}
	}

	/**
	 * Gives same result as Math.toRadians for some particular values like 90.0,
	 * 180.0 or 360.0, but is faster (no division).
	 * 
	 * @param angdeg Angle value in degrees.
	 * @return Angle value in radians.
	 */
	public static double toRadians(double angdeg) {
		if (USE_JDK_MATH) {
			return Math.toRadians(angdeg);
		}
		return angdeg * (Math.PI / 180);
	}

	/**
	 * Gives same result as Math.toDegrees for some particular values like
	 * Math.PI/2, Math.PI or 2*Math.PI, but is faster (no division).
	 * 
	 * @param angrad Angle value in radians.
	 * @return Angle value in degrees.
	 */
	public static double toDegrees(double angrad) {
		if (USE_JDK_MATH) {
			return Math.toDegrees(angrad);
		}
		return angrad * (180 / Math.PI);
	}

	/**
	 * @param sign    Sign of the angle: true for positive, false for negative.
	 * @param degrees Degrees, in [0,180].
	 * @param minutes Minutes, in [0,59].
	 * @param seconds Seconds, in [0.0,60.0[.
	 * @return Angle in radians.
	 */
	public static double toRadians(boolean sign, int degrees, int minutes, double seconds) {
		return toRadians(toDegrees(sign, degrees, minutes, seconds));
	}

	/**
	 * @param sign    Sign of the angle: true for positive, false for negative.
	 * @param degrees Degrees, in [0,180].
	 * @param minutes Minutes, in [0,59].
	 * @param seconds Seconds, in [0.0,60.0[.
	 * @return Angle in degrees.
	 */
	public static double toDegrees(boolean sign, int degrees, int minutes, double seconds) {
		double signFactor = sign ? 1.0 : -1.0;
		return signFactor * (degrees + (1.0 / 60) * (minutes + (1.0 / 60) * seconds));
	}

	/**
	 * @param angrad  Angle in radians.
	 * @param degrees (out) Degrees, in [0,180].
	 * @param minutes (out) Minutes, in [0,59].
	 * @param seconds (out) Seconds, in [0.0,60.0[.
	 * @return true if the resulting angle in [-180deg,180deg] is positive, false if
	 *         it is negative.
	 */
	public static boolean toDMS(double angrad, IntWrapper degrees, IntWrapper minutes, DoubleWrapper seconds) {
		// Computing longitude DMS.
		double tmp = toDegrees(normalizeMinusPiPi(angrad));
		boolean isNeg = (tmp < 0.0);
		if (isNeg) {
			tmp = -tmp;
		}
		degrees.value = (int) tmp;
		tmp = (tmp - degrees.value) * 60.0;
		minutes.value = (int) tmp;
		seconds.value = Math.min((tmp - minutes.value) * 60.0, DOUBLE_BEFORE_60);
		return !isNeg;
	}

	/**
	 * NB: Since 2*Math.PI < 2*PI, a span of 2*Math.PI does not mean full angular
	 * range. ex.: isInClockwiseDomain(0.0, 2*Math.PI, -1e-20) returns false. --->
	 * For full angular range, use a span > 2*Math.PI, like 2*PI_SUP constant of
	 * this class.
	 * 
	 * @param startAngRad An angle, in radians.
	 * @param angSpanRad  An angular span, >= 0.0, in radians.
	 * @param angRad      An angle, in radians.
	 * @return true if angRad is in the clockwise angular domain going from
	 *         startAngRad, over angSpanRad, extremities included, false otherwise.
	 */
	public static boolean isInClockwiseDomain(double startAngRad, double angSpanRad, double angRad) {
		if (Math.abs(angRad) < -TWO_MATH_PI_IN_MINUS_PI_PI) {
			// special case for angular values of small magnitude
			if (angSpanRad <= 2 * Math.PI) {
				if (angSpanRad < 0.0) {
					// empty domain
					return false;
				}
				// angSpanRad is in [0,2*PI]
				startAngRad = normalizeMinusPiPi(startAngRad);
				double endAngRad = normalizeMinusPiPi(startAngRad + angSpanRad);
				if (startAngRad <= endAngRad) {
					return (angRad >= startAngRad) && (angRad <= endAngRad);
				} else {
					return (angRad >= startAngRad) || (angRad <= endAngRad);
				}
			} else { // angSpanRad > 2*Math.PI, or is NaN
				return (angSpanRad == angSpanRad);
			}
		} else {
			// general case
			return (normalizeZeroTwoPi(angRad - startAngRad) <= angSpanRad);
		}
	}

	/*
	 * hyperbolic trigonometry
	 */

	/**
	 * Some properties of sinh(x) = (exp(x)-exp(-x))/2: 1) defined on
	 * ]-Infinity,+Infinity[ 2) result in ]-Infinity,+Infinity[ 3) sinh(x) =
	 * -sinh(-x) (implies sinh(0) = 0) 4) sinh(epsilon) ~= epsilon 5)
	 * lim(sinh(x),x->+Infinity) = +Infinity (y increasing exponentially faster than
	 * x) 6) reaches +Infinity (double overflow) for x >= 710.475860073944, i.e. a
	 * bit further than exp(x)
	 * 
	 * @param value A double value.
	 * @return Value hyperbolic sine.
	 */
	public static double sinh(double value) {
		if (USE_JDK_MATH) {
			return Math.sinh(value);
		}
		// sinh(x) = (exp(x)-exp(-x))/2
		double h;
		if (value < 0.0) {
			value = -value;
			h = -0.5;
		} else {
			h = 0.5;
		}
		if (value < 22.0) {
			if (value < TWO_POW_N28) {
				return (h < 0.0) ? -value : value;
			} else {
				// sinh(x)
				// = (exp(x)-exp(-x))/2
				// = (exp(x)-1/exp(x))/2
				// = (expm1(x) + 1 - 1/(expm1(x)+1))/2
				// = (expm1(x) + (expm1(x)+1)/(expm1(x)+1) - 1/(expm1(x)+1))/2
				// = (expm1(x) + expm1(x)/(expm1(x)+1))/2
				double t = expm1(value);
				// Might be more accurate, if value < 1: return h*((t+t)-t*t/(t+1.0)).
				return h * (t + t / (t + 1.0));
			}
		} else if (value < LOG_DOUBLE_MAX_VALUE) {
			return h * exp(value);
		} else {
			double t = exp(value * 0.5);
			return (h * t) * t;
		}
	}

	/**
	 * Some properties of cosh(x) = (exp(x)+exp(-x))/2: 1) defined on
	 * ]-Infinity,+Infinity[ 2) result in [1,+Infinity[ 3) cosh(0) = 1 4) cosh(x) =
	 * cosh(-x) 5) lim(cosh(x),x->+Infinity) = +Infinity (y increasing exponentially
	 * faster than x) 6) reaches +Infinity (double overflow) for x >=
	 * 710.475860073944, i.e. a bit further than exp(x)
	 * 
	 * @param value A double value.
	 * @return Value hyperbolic cosine.
	 */
	public static double cosh(double value) {
		if (USE_JDK_MATH) {
			return Math.cosh(value);
		}
		// cosh(x) = (exp(x)+exp(-x))/2
		if (value < 0.0) {
			value = -value;
		}
		if (value < LOG_TWO_POW_27) {
			if (value < TWO_POW_N27) {
				// cosh(x)
				// = (exp(x)+exp(-x))/2
				// = ((1+x+x^2/2!+...) + (1-x+x^2/2!-...))/2
				// = 1+x^2/2!+x^4/4!+...
				// For value of x small in magnitude, the sum of the terms does not add to 1.
				return 1;
			} else {
				// cosh(x)
				// = (exp(x)+exp(-x))/2
				// = (exp(x)+1/exp(x))/2
				double t = exp(value);
				return 0.5 * (t + 1 / t);
			}
		} else if (value < LOG_DOUBLE_MAX_VALUE) {
			return 0.5 * exp(value);
		} else {
			double t = exp(value * 0.5);
			return (0.5 * t) * t;
		}
	}

	/**
	 * Much more accurate than cosh(value)-1, for arguments (and results) close to
	 * zero.
	 * 
	 * coshm1(-0.0) = -0.0, for homogeneity with acosh1p(-0.0) = -0.0.
	 * 
	 * @param value A double value.
	 * @return Value hyperbolic cosine, minus 1.
	 */
	public static double coshm1(double value) {
		// cosh(x)-1 = (exp(x)+exp(-x))/2 - 1
		if (value < 0.0) {
			value = -value;
		}
		if (value < LOG_TWO_POW_27) {
			if (value < TWO_POW_N27) {
				if (value == 0.0) {
					// +-0.0
					return value;
				}
				// Using (expm1(x)+expm1(-x))/2
				// is not accurate for tiny values,
				// for expm1 results are of higher
				// magnitude than the result and
				// of different signs, such as their
				// sum is not accurate.
				// cosh(x) - 1
				// = (exp(x)+exp(-x))/2 - 1
				// = ((1+x+x^2/2!+...) + (1-x+x^2/2!-...))/2 - 1
				// = x^2/2!+x^4/4!+...
				// ~= x^2 * (1/2 + x^2 * 1/24)
				// = x^2 * 0.5 (since x < 2^-27)
				return 0.5 * value * value;
			} else {
				// cosh(x) - 1
				// = (exp(x)+exp(-x))/2 - 1
				// = (exp(x)-1+exp(-x)-1)/2
				// = (expm1(x)+expm1(-x))/2
				return 0.5 * (expm1(value) + expm1(-value));
			}
		} else if (value < LOG_DOUBLE_MAX_VALUE) {
			return 0.5 * exp(value) - 1.0;
		} else {
			// No need to subtract 1 from result.
			double t = exp(value * 0.5);
			return (0.5 * t) * t;
		}
	}

	/**
	 * Computes hyperbolic sine and hyperbolic cosine together.
	 * 
	 * @param value   A double value.
	 * @param hcosine (out) Value hyperbolic cosine.
	 * @return Value hyperbolic sine.
	 */
	public static double sinhAndCosh(double value, DoubleWrapper hcosine) {
		if (USE_JDK_MATH) {
			hcosine.value = Math.cosh(value);
			return Math.sinh(value);
		}
		// Mixup of sinh and cosh treatments: if you modify them,
		// you might want to also modify this.
		double h;
		if (value < 0.0) {
			value = -value;
			h = -0.5;
		} else {
			h = 0.5;
		}
		final double hsine;
		// LOG_TWO_POW_27 = 18.714973875118524
		if (value < LOG_TWO_POW_27) { // test from cosh
			// sinh
			if (value < TWO_POW_N28) {
				hsine = (h < 0.0) ? -value : value;
			} else {
				double t = expm1(value);
				hsine = h * (t + t / (t + 1.0));
			}
			// cosh
			if (value < TWO_POW_N27) {
				hcosine.value = 1;
			} else {
				double t = exp(value);
				hcosine.value = 0.5 * (t + 1 / t);
			}
		} else if (value < 22.0) { // test from sinh
			// Here, value is in [18.714973875118524,22.0[.
			double t = expm1(value);
			hsine = h * (t + t / (t + 1.0));
			hcosine.value = 0.5 * (t + 1.0);
		} else {
			if (value < LOG_DOUBLE_MAX_VALUE) {
				hsine = h * exp(value);
			} else {
				double t = exp(value * 0.5);
				hsine = (h * t) * t;
			}
			hcosine.value = Math.abs(hsine);
		}
		return hsine;
	}

	/**
	 * Some properties of tanh(x) = sinh(x)/cosh(x) = (exp(2*x)-1)/(exp(2*x)+1): 1)
	 * defined on ]-Infinity,+Infinity[ 2) result in ]-1,1[ 3) tanh(x) = -tanh(-x)
	 * (implies tanh(0) = 0) 4) tanh(epsilon) ~= epsilon 5)
	 * lim(tanh(x),x->+Infinity) = 1 6) reaches 1 (double loss of precision) for x =
	 * 19.061547465398498
	 * 
	 * @param value A double value.
	 * @return Value hyperbolic tangent.
	 */
	public static double tanh(double value) {
		if (USE_JDK_MATH) {
			return Math.tanh(value);
		}
		// tanh(x) = sinh(x)/cosh(x)
		// = (exp(x)-exp(-x))/(exp(x)+exp(-x))
		// = (exp(2*x)-1)/(exp(2*x)+1)
		boolean negateResult = false;
		if (value < 0.0) {
			value = -value;
			negateResult = true;
		}
		double z;
		if (value < TANH_1_THRESHOLD) {
			if (value < TWO_POW_N55) {
				return negateResult ? -value * (1.0 - value) : value * (1.0 + value);
			} else if (value >= 1) {
				z = 1.0 - 2.0 / (expm1(value + value) + 2.0);
			} else {
				double t = expm1(-(value + value));
				z = -t / (t + 2.0);
			}
		} else {
			z = (value != value) ? Double.NaN : 1.0;
		}
		return negateResult ? -z : z;
	}

	/**
	 * Some properties of asinh(x) = log(x + sqrt(x^2 + 1)) 1) defined on
	 * ]-Infinity,+Infinity[ 2) result in ]-Infinity,+Infinity[ 3) asinh(x) =
	 * -asinh(-x) (implies asinh(0) = 0) 4) asinh(epsilon) ~= epsilon 5)
	 * lim(asinh(x),x->+Infinity) = +Infinity (y increasing logarithmically slower
	 * than x)
	 * 
	 * @param value A double value.
	 * @return Value hyperbolic arcsine.
	 */
	public static double asinh(double value) {
		// asinh(x) = log(x + sqrt(x^2 + 1))
		boolean negateResult = false;
		if (value < 0.0) {
			value = -value;
			negateResult = true;
		}
		double result;
		// (about) smallest possible for
		// non-log1p case to be accurate.
		if (value < ASINH_LOG1P_THRESHOLD) {
			// Around this range, FDLIBM uses
			// log1p(value+value*value/(1+sqrt(value*value+1))),
			// but it's slower, so we don't use it.
			/*
			 * If x is close to zero, log argument is close to 1, so to avoid precision loss
			 * we use log1p(double), with (1+x)^p = 1 + p * x + (p*(p-1))/2! * x^2 +
			 * (p*(p-1)*(p-2))/3! * x^3 + ... (1+x)^p = 1 + p * x * (1 + (p-1)/2 * x * (1 +
			 * (p-2)/3 * x + ...) (1+x)^0.5 = 1 + 0.5 * x * (1 + (0.5-1)/2 * x * (1 +
			 * (0.5-2)/3 * x + ...) (1+x^2)^0.5 = 1 + 0.5 * x^2 * (1 + (0.5-1)/2 * x^2 * (1
			 * + (0.5-2)/3 * x^2 + ...) x + (1+x^2)^0.5 = 1 + x * (1 + 0.5 * x * (1 +
			 * (0.5-1)/2 * x^2 * (1 + (0.5-2)/3 * x^2 + ...)) so asinh(x) = log1p(x * (1 +
			 * 0.5 * x * (1 + (0.5-1)/2 * x^2 * (1 + (0.5-2)/3 * x^2 + ...)))
			 */
			final double x = value;
			final double x2 = x * x;
			// Enough terms for good accuracy,
			// given our threshold.
			final double argLog1p = (x
					* (1 + 0.5 * x * (1 + (0.5 - 1) / 2 * x2 * (1 + (0.5 - 2) / 3 * x2 * (1 + (0.5 - 3) / 4 * x2 * (1 + (0.5 - 4) / 5 * x2))))));
			result = log1p(argLog1p);
		} else if (value < ASINH_ACOSH_SQRT_ELISION_THRESHOLD) {
			// Around this range, FDLIBM uses
			// log(2*value+1/(value+sqrt(value*value+1))),
			// but it involves an additional division
			// so we don't use it.
			result = log(value + sqrt(value * value + 1.0));
		} else {
			// log(2*value) would overflow for value > Double.MAX_VALUE/2,
			// so we compute otherwise.
			result = LOG_2 + log(value);
		}
		return negateResult ? -result : result;
	}

	/**
	 * Some properties of acosh(x) = log(x + sqrt(x^2 - 1)): 1) defined on
	 * [1,+Infinity[ 2) result in ]0,+Infinity[ (by convention, since cosh(x) =
	 * cosh(-x)) 3) acosh(1) = 0 4) acosh(1+epsilon) ~= log(1 + sqrt(2*epsilon)) ~=
	 * sqrt(2*epsilon) 5) lim(acosh(x),x->+Infinity) = +Infinity (y increasing
	 * logarithmically slower than x)
	 * 
	 * @param value A double value.
	 * @return Value hyperbolic arccosine.
	 */
	public static double acosh(double value) {
		if (!(value > 1.0)) {
			// NaN, or value <= 1
			if (ANTI_JIT_OPTIM_CRASH_ON_NAN) {
				return (value < 1.0) ? Double.NaN : value - 1.0;
			} else {
				return (value == 1.0) ? 0.0 : Double.NaN;
			}
		}
		double result;
		if (value < ASINH_ACOSH_SQRT_ELISION_THRESHOLD) {
			// Around this range, FDLIBM uses
			// log(2*value-1/(value+sqrt(value*value-1))),
			// but it involves an additional division
			// so we don't use it.
			result = log(value + sqrt(value * value - 1.0));
		} else {
			// log(2*value) would overflow for value > Double.MAX_VALUE/2,
			// so we compute otherwise.
			result = LOG_2 + log(value);
		}
		return result;
	}

	/**
	 * Much more accurate than acosh(1+value), for arguments (and results) close to
	 * zero.
	 * 
	 * acosh1p(-0.0) = -0.0, for homogeneity with sqrt(-0.0) = -0.0, which looks
	 * about the same near 0.
	 * 
	 * @param value A double value.
	 * @return Hyperbolic arccosine of (1+value).
	 */
	public static double acosh1p(double value) {
		if (!(value > 0.0)) {
			// NaN, or value <= 0.
			// If value is -0.0, returning -0.0.
			if (ANTI_JIT_OPTIM_CRASH_ON_NAN) {
				return (value < 0.0) ? Double.NaN : value;
			} else {
				return (value == 0.0) ? value : Double.NaN;
			}
		}
		double result;
		if (value < (ASINH_ACOSH_SQRT_ELISION_THRESHOLD - 1)) {
			// acosh(1+x)
			// = log((1+x) + sqrt((1+x)^2 - 1))
			// = log(1 + x + sqrt(1 + 2*x + x^2 - 1))
			// = log1p(x + sqrt(2*x + x^2))
			// = log1p(x + sqrt(x * (2 + x))
			result = log1p(value + sqrt(value * (2 + value)));
		} else {
			result = LOG_2 + log(1 + value);
		}
		return result;
	}

	/**
	 * Some properties of atanh(x) = log((1+x)/(1-x))/2: 1) defined on ]-1,1[ 2)
	 * result in ]-Infinity,+Infinity[ 3) atanh(-1) = -Infinity (by continuity) 4)
	 * atanh(1) = +Infinity (by continuity) 5) atanh(epsilon) ~= epsilon 6)
	 * lim(atanh(x),x->1) = +Infinity
	 * 
	 * @param value A double value.
	 * @return Value hyperbolic arctangent.
	 */
	public static double atanh(double value) {
		boolean negateResult = false;
		if (value < 0.0) {
			value = -value;
			negateResult = true;
		}
		double result;
		if (!(value < 1.0)) {
			// NaN, or value >= 1
			if (ANTI_JIT_OPTIM_CRASH_ON_NAN) {
				result = (value > 1.0) ? Double.NaN : Double.POSITIVE_INFINITY + value;
			} else {
				result = (value == 1.0) ? Double.POSITIVE_INFINITY : Double.NaN;
			}
		} else {
			// For value < 0.5, FDLIBM uses
			// 0.5 * log1p((value+value) + (value+value)*value/(1-value)),
			// instead, but this is good enough for us.
			// atanh(x)
			// = log((1+x)/(1-x))/2
			// = log((1-x+2x)/(1-x))/2
			// = log1p(2x/(1-x))/2
			result = 0.5 * log1p((value + value) / (1.0 - value));
		}
		return negateResult ? -result : result;
	}

	/*
	 * exponentials
	 */

	/**
	 * @param value A double value.
	 * @return e^value.
	 */
	public static double exp(double value) {
		if (USE_JDK_MATH) {
			return Math.exp(value);
		}
		// exp(x) = exp([x])*exp(y)
		// with [x] the integer part of x, and y = x-[x]
		// ===>
		// We find an approximation of y, called z.
		// ===>
		// exp(x) = exp([x])*(exp(z)*exp(epsilon))
		// with epsilon = y - z
		// ===>
		// We have exp([x]) and exp(z) pre-computed in tables, we "just" have to compute
		// exp(epsilon).
		//
		// We use the same indexing (cast to int) to compute x integer part and the
		// table index corresponding to z, to avoid two int casts.
		// Also, to optimize index multiplication and division, we use powers of two,
		// so that we can do it with bits shifts.

		if (value > EXP_OVERFLOW_LIMIT) {
			return Double.POSITIVE_INFINITY;
		} else if (!(value >= EXP_UNDERFLOW_LIMIT)) {
			return (value != value) ? Double.NaN : 0.0;
		}

		final int indexes = (int) (value * EXP_LO_INDEXING);

		final int valueInt;
		if (indexes >= 0) {
			valueInt = (indexes >> EXP_LO_INDEXING_DIV_SHIFT);
		} else {
			valueInt = -((-indexes) >> EXP_LO_INDEXING_DIV_SHIFT);
		}
		final double hiTerm = MyTExp.expHiTab[valueInt - (int) EXP_UNDERFLOW_LIMIT];

		final int zIndex = indexes - (valueInt << EXP_LO_INDEXING_DIV_SHIFT);
		final double y = (value - valueInt);
		final double z = zIndex * (1.0 / EXP_LO_INDEXING);
		final double eps = y - z;
		final double expZ = MyTExp.expLoPosTab[zIndex + EXP_LO_TAB_MID_INDEX];
		final double expEps = (1 + eps * (1 + eps * (1.0 / 2 + eps * (1.0 / 6 + eps * (1.0 / 24)))));
		final double loTerm = expZ * expEps;

		return hiTerm * loTerm;
	}

	/**
	 * Quick exp, with a max relative error of about 2.94e-2 for |value| < 700.0 or
	 * so, and no accuracy at all outside this range. Derived from a note by Nicol
	 * N. Schraudolph, IDSIA, 1998.
	 * 
	 * @param value A double value.
	 * @return e^value.
	 */
	public static double expQuick(double value) {
		if (USE_JDK_MATH) {
			return Math.exp(value);
		}
		/*
		 * Cast of double values, even in long range, into long, is slower than from
		 * double to int for values in int range, and then from int to long. For that
		 * reason, we only work with integer values in int range (corresponding to the
		 * 32 first bits of the long, containing sign, exponent, and highest significant
		 * bits of double's mantissa), and cast twice.
		 * 
		 * Constants determined empirically, using a random-based metaheuristic. Should
		 * be possible to find better ones.
		 */
		return Double.longBitsToDouble(((long) (int) (1512775.3952 * value + 1.0726481222E9)) << 32);
	}

	/**
	 * Much more accurate than exp(value)-1, for arguments (and results) close to
	 * zero.
	 * 
	 * @param value A double value.
	 * @return e^value-1.
	 */
	public static double expm1(double value) {
		if (USE_JDK_MATH) {
			return Math.expm1(value);
		}
		// If value is far from zero, we use exp(value)-1.
		//
		// If value is close to zero, we use the following formula:
		// exp(value)-1
		// = exp(valueApprox)*exp(epsilon)-1
		// = exp(valueApprox)*(exp(epsilon)-exp(-valueApprox))
		// = exp(valueApprox)*(1+epsilon+epsilon^2/2!+...-exp(-valueApprox))
		// = exp(valueApprox)*((1-exp(-valueApprox))+epsilon+epsilon^2/2!+...)
		// exp(valueApprox) and exp(-valueApprox) being stored in tables.

		if (Math.abs(value) < EXP_LO_DISTANCE_TO_ZERO) {
			// Taking int part instead of rounding, which takes too long.
			int i = (int) (value * EXP_LO_INDEXING);
			double delta = value - i * (1.0 / EXP_LO_INDEXING);
			return MyTExp.expLoPosTab[i + EXP_LO_TAB_MID_INDEX] * (MyTExp.expLoNegTab[i + EXP_LO_TAB_MID_INDEX]
					+ delta * (1 + delta * (1.0 / 2 + delta * (1.0 / 6 + delta * (1.0 / 24 + delta * (1.0 / 120))))));
		} else {
			return exp(value) - 1;
		}
	}

	/*
	 * logarithms
	 */

	/**
	 * @param value A double value.
	 * @return Value logarithm (base e).
	 */
	public static double log(double value) {
		if (USE_JDK_MATH || (!USE_REDEFINED_LOG)) {
			return Math.log(value);
		}
		if (value > 0.0) {
			if (value == Double.POSITIVE_INFINITY) {
				return Double.POSITIVE_INFINITY;
			}

			// For normal values not close to 1.0, we use the following formula:
			// log(value)
			// = log(2^exponent*1.mantissa)
			// = log(2^exponent) + log(1.mantissa)
			// = exponent * log(2) + log(1.mantissa)
			// = exponent * log(2) + log(1.mantissaApprox) +
			// log(1.mantissa/1.mantissaApprox)
			// = exponent * log(2) + log(1.mantissaApprox) + log(1+epsilon)
			// = exponent * log(2) + log(1.mantissaApprox) +
			// epsilon-epsilon^2/2+epsilon^3/3-epsilon^4/4+...
			// with:
			// 1.mantissaApprox <= 1.mantissa,
			// log(1.mantissaApprox) in table,
			// epsilon = (1.mantissa/1.mantissaApprox)-1
			//
			// To avoid bad relative error for small results,
			// values close to 1.0 are treated aside, with the formula:
			// log(x) = z*(2+z^2*((2.0/3)+z^2*((2.0/5))+z^2*((2.0/7))+...)))
			// with z=(x-1)/(x+1)

			double h;
			if (value > 0.95) {
				if (value < 1.14) {
					double z = (value - 1.0) / (value + 1.0);
					double z2 = z * z;
					return z * (2 + z2 * ((2.0 / 3) + z2 * ((2.0 / 5) + z2 * ((2.0 / 7) + z2 * ((2.0 / 9) + z2 * ((2.0 / 11)))))));
				}
				h = 0.0;
			} else if (value < DOUBLE_MIN_NORMAL) {
				// Ensuring value is normal.
				value *= TWO_POW_52;
				// log(x*2^52)
				// = log(x)-ln(2^52)
				// = log(x)-52*ln(2)
				h = -52 * LOG_2;
			} else {
				h = 0.0;
			}

			int valueBitsHi = (int) (Double.doubleToRawLongBits(value) >> 32);
			int valueExp = (valueBitsHi >> 20) - MAX_DOUBLE_EXPONENT;
			// Getting the first LOG_BITS bits of the mantissa.
			int xIndex = ((valueBitsHi << 12) >>> (32 - LOG_BITS));

			// 1.mantissa/1.mantissaApprox - 1
			double z = (value * twoPowNormalOrSubnormal(-valueExp)) * MyTLog.logXInvTab[xIndex] - 1;

			z *= (1 - z * ((1.0 / 2) - z * ((1.0 / 3))));

			return h + valueExp * LOG_2 + (MyTLog.logXLogTab[xIndex] + z);

		} else if (value == 0.0) {
			return Double.NEGATIVE_INFINITY;
		} else { // value < 0.0, or value is NaN
			return Double.NaN;
		}
	}

	/**
	 * Quick log, with a max relative error of about 1.9e-3 for values in
	 * ]Double.MIN_NORMAL,+Infinity[, and worse accuracy outside this range.
	 * 
	 * @param value A double value, in ]0,+Infinity[ (strictly positive and finite).
	 * @return Value logarithm (base e).
	 */
	public static double logQuick(double value) {
		if (USE_JDK_MATH) {
			return Math.log(value);
		}
		/*
		 * Inverse of Schraudolph's method for exp, is very inaccurate near 1, and not
		 * that fast (even using floats), especially with added if's to deal with values
		 * near 1, so we don't use it, and use a simplified version of our log's
		 * redefined algorithm.
		 */

		// Simplified version of log's redefined algorithm:
		// log(value) ~= exponent * log(2) + log(1.mantissaApprox)

		double h;
		if (value > 0.87) {
			if (value < 1.16) {
				return 2.0 * (value - 1.0) / (value + 1.0);
			}
			h = 0.0;
		} else if (value < DOUBLE_MIN_NORMAL) {
			value *= TWO_POW_52;
			h = -52 * LOG_2;
		} else {
			h = 0.0;
		}

		int valueBitsHi = (int) (Double.doubleToRawLongBits(value) >> 32);
		int valueExp = (valueBitsHi >> 20) - MAX_DOUBLE_EXPONENT;
		int xIndex = ((valueBitsHi << 12) >>> (32 - LOG_BITS));

		return h + valueExp * LOG_2 + MyTLog.logXLogTab[xIndex];
	}

	/**
	 * @param value A double value.
	 * @return Value logarithm (base 10).
	 */
	public static double log10(double value) {
		if (USE_JDK_MATH || (!USE_REDEFINED_LOG)) {
			return Math.log10(value);
		}
		// INV_LOG_10 is < 1, but there is no risk of log(double)
		// overflow (positive or negative) while the end result shouldn't,
		// since log(Double.MIN_VALUE) and log(Double.MAX_VALUE) have
		// magnitudes of just a few hundreds.
		return log(value) * INV_LOG_10;
	}

	/**
	 * Much more accurate than log(1+value), for arguments (and results) close to
	 * zero.
	 * 
	 * @param value A double value.
	 * @return Logarithm (base e) of (1+value).
	 */
	public static double log1p(double value) {
		if (USE_JDK_MATH) {
			return Math.log1p(value);
		}
		if (false) {
			// This also works. Simpler but a bit slower.
			if (value == Double.POSITIVE_INFINITY) {
				return Double.POSITIVE_INFINITY;
			}
			double valuePlusOne = 1 + value;
			if (valuePlusOne == 1.0) {
				return value;
			} else {
				return log(valuePlusOne) * (value / (valuePlusOne - 1.0));
			}
		}
		if (value > -1.0) {
			if (value == Double.POSITIVE_INFINITY) {
				return Double.POSITIVE_INFINITY;
			}

			// ln'(x) = 1/x
			// so
			// log(x+epsilon) ~= log(x) + epsilon/x
			//
			// Let u be 1+value rounded:
			// 1+value = u+epsilon
			//
			// log(1+value)
			// = log(u+epsilon)
			// ~= log(u) + epsilon/value
			// We compute log(u) as done in log(double), and then add the corrective term.

			double valuePlusOne = 1.0 + value;
			if (valuePlusOne == 1.0) {
				return value;
			} else if (Math.abs(value) < 0.15) {
				double z = value / (value + 2.0);
				double z2 = z * z;
				return z * (2 + z2 * ((2.0 / 3) + z2 * ((2.0 / 5) + z2 * ((2.0 / 7) + z2 * ((2.0 / 9) + z2 * ((2.0 / 11)))))));
			}

			int valuePlusOneBitsHi = (int) (Double.doubleToRawLongBits(valuePlusOne) >> 32) & 0x7FFFFFFF;
			int valuePlusOneExp = (valuePlusOneBitsHi >> 20) - MAX_DOUBLE_EXPONENT;
			// Getting the first LOG_BITS bits of the mantissa.
			int xIndex = ((valuePlusOneBitsHi << 12) >>> (32 - LOG_BITS));

			// 1.mantissa/1.mantissaApprox - 1
			double z = (valuePlusOne * twoPowNormalOrSubnormal(-valuePlusOneExp)) * MyTLog.logXInvTab[xIndex] - 1;

			z *= (1 - z * ((1.0 / 2) - z * (1.0 / 3)));

			// Adding epsilon/valuePlusOne to z,
			// with
			// epsilon = value - (valuePlusOne-1)
			// (valuePlusOne + epsilon ~= 1+value (not rounded))

			return valuePlusOneExp * LOG_2 + MyTLog.logXLogTab[xIndex] + (z + (value - (valuePlusOne - 1)) / valuePlusOne);
		} else if (value == -1.0) {
			return Double.NEGATIVE_INFINITY;
		} else { // value < -1.0, or value is NaN
			return Double.NaN;
		}
	}

	/*
	 * powers
	 */

	/**
	 * 1e-13ish accuracy or better on whole double range.
	 * 
	 * @param value A double value.
	 * @param power A power.
	 * @return value^power.
	 */
	public static double pow(double value, double power) {
		if (USE_JDK_MATH) {
			return Math.pow(value, power);
		}
		if (power == 0.0) {
			return 1.0;
		} else if (power == 1.0) {
			return value;
		}
		if (value <= 0.0) {
			// powerInfo: 0 if not integer, 1 if even integer, -1 if odd integer
			int powerInfo;
			if (Math.abs(power) >= (TWO_POW_52 * 2)) {
				// The binary digit just before comma is outside mantissa,
				// thus it is always 0: power is an even integer.
				powerInfo = 1;
			} else {
				// If power's magnitude permits, we cast into int instead of into long,
				// as it is faster.
				if (Math.abs(power) <= (double) Integer.MAX_VALUE) {
					int powerAsInt = (int) power;
					if (power == (double) powerAsInt) {
						powerInfo = ((powerAsInt & 1) == 0) ? 1 : -1;
					} else { // power is not an integer (and not NaN, due to test against Integer.MAX_VALUE)
						powerInfo = 0;
					}
				} else {
					long powerAsLong = (long) power;
					if (power == (double) powerAsLong) {
						powerInfo = ((powerAsLong & 1) == 0) ? 1 : -1;
					} else { // power is not an integer, or is NaN
						if (power != power) {
							return Double.NaN;
						}
						powerInfo = 0;
					}
				}
			}

			if (value == 0.0) {
				if (power < 0.0) {
					return (powerInfo < 0) ? 1 / value : Double.POSITIVE_INFINITY;
				} else { // power > 0.0 (0 and NaN cases already treated)
					return (powerInfo < 0) ? value : 0.0;
				}
			} else { // value < 0.0
				if (value == Double.NEGATIVE_INFINITY) {
					if (powerInfo < 0) { // power odd integer
						return (power < 0.0) ? -0.0 : Double.NEGATIVE_INFINITY;
					} else { // power even integer, or not an integer
						return (power < 0.0) ? 0.0 : Double.POSITIVE_INFINITY;
					}
				} else {
					return (powerInfo == 0) ? Double.NaN : powerInfo * exp(power * log(-value));
				}
			}
		} else { // value > 0.0, or value is NaN
			return exp(power * log(value));
		}
	}

	/**
	 * Quick pow, with a max relative error of about 1e-2 for value >=
	 * Double.MIN_NORMAL and 1e-10 < |value^power| < 1e10, of about 6e-2 for value
	 * >= Double.MIN_NORMAL and 1e-40 < |value^power| < 1e40, and worse accuracy
	 * otherwise.
	 * 
	 * @param value A double value, in ]0,+Infinity[ (strictly positive and finite).
	 * @param power A double value.
	 * @return value^power.
	 */
	public static double powQuick(double value, double power) {
		if (USE_JDK_MATH) {
			return Math.pow(value, power);
		}
		return exp(power * logQuick(value));
	}

	/**
	 * This treatment is somehow accurate for low values of |power|, and for
	 * |power*getExponent(value)| < 1023 or so (to stay away from double extreme
	 * magnitudes (large and small)).
	 * 
	 * @param value A double value.
	 * @param power A power.
	 * @return value^power.
	 */
	public static double powFast(double value, int power) {
		if (USE_JDK_MATH) {
			return Math.pow(value, power);
		}
		if (power < 3) {
			if (power < 0) {
				// Opposite of Integer.MIN_VALUE does not exist as int.
				if (power == Integer.MIN_VALUE) {
					// Integer.MAX_VALUE = -(power+1)
					return 1.0 / (powFast(value, Integer.MAX_VALUE) * value);
				} else {
					return 1.0 / powFast(value, -power);
				}
			} else {
				// Here, power is in [0,2].
				if (power == 2) { // Most common case first.
					return value * value;
				} else if (power == 0) {
					return 1.0;
				} else { // power == 1
					return value;
				}
			}
		} else { // power >= 4
			double oddRemains = 1.0;
			// If power <= 5, faster to finish outside the loop.
			while (power > 5) {
				// Test if power is odd.
				if ((power & 1) != 0) {
					oddRemains *= value;
				}
				value *= value;
				power >>= 1; // power = power / 2
			}
			// Here, power is in [3,5].
			if (power == 3) {
				return oddRemains * value * value * value;
			} else { // power in [4,5].
				double v2 = value * value;
				if (power == 4) {
					return oddRemains * v2 * v2;
				} else { // power == 5
					return oddRemains * v2 * v2 * value;
				}
			}
		}
	}

	/**
	 * @param value A float value.
	 * @return value*value.
	 */
	public static float pow2(float value) {
		return value * value;
	}

	/**
	 * @param value A double value.
	 * @return value*value.
	 */
	public static double pow2(double value) {
		return value * value;
	}

	/**
	 * @param value A float value.
	 * @return value*value*value.
	 */
	public static float pow3(float value) {
		return value * value * value;
	}

	/**
	 * @param value A double value.
	 * @return value*value*value.
	 */
	public static double pow3(double value) {
		return value * value * value;
	}

	/*
	 * roots
	 */

	/**
	 * @param value A double value.
	 * @return Value square root.
	 */
	public static double sqrt(double value) {
		if (USE_JDK_MATH || (!USE_REDEFINED_SQRT)) {
			return Math.sqrt(value);
		}
		// See cbrt for comments, sqrt uses the same ideas.

		if (!(value > 0.0)) { // value <= 0.0, or value is NaN
			if (ANTI_JIT_OPTIM_CRASH_ON_NAN) {
				return (value < 0.0) ? Double.NaN : value;
			} else {
				return (value == 0.0) ? value : Double.NaN;
			}
		} else if (value == Double.POSITIVE_INFINITY) {
			return Double.POSITIVE_INFINITY;
		}

		double h;
		if (value < DOUBLE_MIN_NORMAL) {
			value *= TWO_POW_52;
			h = 2 * TWO_POW_N26;
		} else {
			h = 2.0;
		}

		int valueBitsHi = (int) (Double.doubleToRawLongBits(value) >> 32);
		int valueExponentIndex = (valueBitsHi >> 20) + (-MAX_DOUBLE_EXPONENT - MIN_DOUBLE_EXPONENT);
		int xIndex = ((valueBitsHi << 12) >>> (32 - SQRT_LO_BITS));

		double result = MyTSqrt.sqrtXSqrtHiTab[valueExponentIndex] * MyTSqrt.sqrtXSqrtLoTab[xIndex];
		double slope = MyTSqrt.sqrtSlopeHiTab[valueExponentIndex] * MyTSqrt.sqrtSlopeLoTab[xIndex];
		value *= 0.25;

		result += (value - result * result) * slope;
		result += (value - result * result) * slope;
		return h * (result + (value - result * result) * slope);
	}

	/**
	 * Quick sqrt, with with a max relative error of about 3.41e-2 for values in
	 * [Double.MIN_NORMAL,Double.MAX_VALUE], and worse accuracy outside this range.
	 * 
	 * @param value A double value.
	 * @return Value square root.
	 */
	public static double sqrtQuick(double value) {
		if (USE_JDK_MATH) {
			return Math.sqrt(value);
		}
		final long bits = Double.doubleToRawLongBits(value);
		/*
		 * Constant determined empirically, using a random-based metaheuristic. Should
		 * be possible to find a better one.
		 */
		return Double.longBitsToDouble((bits + 4606859074900000000L) >>> 1);
	}

	/**
	 * Quick inverse of square root, with a max relative error of about 3.44e-2 for
	 * values in [Double.MIN_NORMAL,Double.MAX_VALUE], and worse accuracy outside
	 * this range.
	 * 
	 * This implementation uses zero step of Newton's method. Here are the max
	 * relative errors on [Double.MIN_NORMAL,Double.MAX_VALUE] depending on number
	 * of steps, if you want to copy-paste this code and use your own number: n=0:
	 * about 3.44e-2 n=1: about 1.75e-3 n=2: about 4.6e-6 n=3: about 3.17e-11 n=4:
	 * about 3.92e-16 n=5: about 3.03e-16
	 * 
	 * @param value A double value.
	 * @return Inverse of value square root.
	 */
	public static double invSqrtQuick(double value) {
		if (USE_JDK_MATH) {
			return 1 / Math.sqrt(value);
		}
		/*
		 * http://en.wikipedia.org/wiki/Fast_inverse_square_root
		 */
		if (false) {
			// With one Newton step (much slower than
			// 1/Math.sqrt(double) if not optimized).
			final double halfInitial = value * 0.5;
			long bits = Double.doubleToRawLongBits(value);
			// If n=0, 6910474759270000000L might be better (3.38e-2 max relative error).
			bits = 0x5FE6EB50C7B537A9L - (bits >> 1);
			value = Double.longBitsToDouble(bits);
			value = value * (1.5 - halfInitial * value * value); // Newton step, can repeat.
			return value;
		} else {
			return Double.longBitsToDouble(0x5FE6EB50C7B537A9L - (Double.doubleToRawLongBits(value) >> 1));
		}
	}

	/**
	 * @param value A double value.
	 * @return Value cubic root.
	 */
	public static double cbrt(double value) {
		if (USE_JDK_MATH) {
			return Math.cbrt(value);
		}
		double h;
		if (value < 0.0) {
			if (value == Double.NEGATIVE_INFINITY) {
				return Double.NEGATIVE_INFINITY;
			}
			value = -value;
			// Making sure value is normal.
			if (value < DOUBLE_MIN_NORMAL) {
				value *= (TWO_POW_52 * TWO_POW_26);
				// h = <result_sign> * <result_multiplicator_to_avoid_overflow> /
				// <cbrt(value_multiplicator_to_avoid_subnormal)>
				h = -2 * TWO_POW_N26;
			} else {
				h = -2.0;
			}
		} else {
			if (!(value < Double.POSITIVE_INFINITY)) { // value is +Infinity, or value is NaN
				return value;
			}
			// Making sure value is normal.
			if (value < DOUBLE_MIN_NORMAL) {
				if (value == 0.0) {
					// cbrt(0.0) = 0.0, cbrt(-0.0) = -0.0
					return value;
				}
				value *= (TWO_POW_52 * TWO_POW_26);
				h = 2 * TWO_POW_N26;
			} else {
				h = 2.0;
			}
		}

		// Normal value is (2^<value exponent> * <a value in [1,2[>).
		// First member cubic root is computed, and multiplied with an approximation
		// of the cubic root of the second member, to end up with a good guess of
		// the result before using Newton's (or Archimedes's) method.
		// To compute the cubic root approximation, we use the formula "cbrt(value) =
		// cbrt(x) * cbrt(value/x)",
		// choosing x as close to value as possible but inferior to it, so that
		// cbrt(value/x) is close to 1
		// (we could iterate on this method, using value/x as new value for each
		// iteration,
		// but finishing with Newton's method is faster).

		// Shift and cast into an int, which overall is faster than working with a long.
		int valueBitsHi = (int) (Double.doubleToRawLongBits(value) >> 32);
		int valueExponentIndex = (valueBitsHi >> 20) + (-MAX_DOUBLE_EXPONENT - MIN_DOUBLE_EXPONENT);
		// Getting the first CBRT_LO_BITS bits of the mantissa.
		int xIndex = ((valueBitsHi << 12) >>> (32 - CBRT_LO_BITS));
		double result = MyTCbrt.cbrtXCbrtHiTab[valueExponentIndex] * MyTCbrt.cbrtXCbrtLoTab[xIndex];
		double slope = MyTCbrt.cbrtSlopeHiTab[valueExponentIndex] * MyTCbrt.cbrtSlopeLoTab[xIndex];

		// Lowering values to avoid overflows when using Newton's method
		// (we will then just have to return twice the result).
		// result^3 = value
		// (result/2)^3 = value/8
		value *= 0.125;
		// No need to divide result here, as division is factorized in result
		// computation tables.
		// result *= 0.5;

		// Newton's method, looking for y = x^(1/p):
		// y(n) = y(n-1) + (x-y(n-1)^p) * slope(y(n-1))
		// y(n) = y(n-1) + (x-y(n-1)^p) * (1/p)*(x(n-1)^(1/p-1))
		// y(n) = y(n-1) + (x-y(n-1)^p) * (1/p)*(x(n-1)^((1-p)/p))
		// with x(n-1)=y(n-1)^p, i.e.:
		// y(n) = y(n-1) + (x-y(n-1)^p) * (1/p)*(y(n-1)^(1-p))
		//
		// For p=3:
		// y(n) = y(n-1) + (x-y(n-1)^3) * (1/(3*y(n-1)^2))

		// To save time, we don't recompute the slope between Newton's method steps,
		// as initial slope is good enough for a few iterations.
		//
		// NB: slope = 1/(3*trueResult*trueResult)
		// As we have result = trueResult/2 (to avoid overflows), we have:
		// slope = 4/(3*result*result)
		// = (4/3)*resultInv*resultInv
		// with newResultInv = 1/newResult
		// = 1/(oldResult+resultDelta)
		// = (oldResultInv)*1/(1+resultDelta/oldResult)
		// = (oldResultInv)*1/(1+resultDelta*oldResultInv)
		// ~= (oldResultInv)*(1-resultDelta*oldResultInv)
		// ===> Successive slopes could be computed without division, if needed,
		// by computing resultInv (instead of slope right away) and retrieving
		// slopes from it.

		result += (value - result * result * result) * slope;
		result += (value - result * result * result) * slope;
		return h * (result + (value - result * result * result) * slope);
	}

	/**
	 * @return sqrt(x^2+y^2) without intermediate overflow or underflow.
	 */
	public static double hypot(double x, double y) {
		if (USE_JDK_MATH) {
			return Math.hypot(x, y);
		}
		x = Math.abs(x);
		y = Math.abs(y);
		// Ensuring x <= y.
		if (y < x) {
			double a = x;
			x = y;
			y = a;
		} else if (!(y >= x)) { // Testing if we have some NaN.
			return hypot_NaN(x, y);
		}

		if (y - x == y) {
			// x too small to subtract from y.
			return y;
		} else {
			double factor;
			if (y > HYPOT_MAX_MAG) {
				// y is too large: scaling down.
				x *= (1 / HYPOT_FACTOR);
				y *= (1 / HYPOT_FACTOR);
				factor = HYPOT_FACTOR;
			} else if (x < (1 / HYPOT_MAX_MAG)) {
				// x is too small: scaling up.
				x *= HYPOT_FACTOR;
				y *= HYPOT_FACTOR;
				factor = (1 / HYPOT_FACTOR);
			} else {
				factor = 1.0;
			}
			return factor * sqrt(x * x + y * y);
		}
	}

	/**
	 * @return sqrt(x^2+y^2+z^2) without intermediate overflow or underflow.
	 */
	public static double hypot(double x, double y, double z) {
		if (USE_JDK_MATH) {
			// No simple JDK equivalent.
		}
		x = Math.abs(x);
		y = Math.abs(y);
		z = Math.abs(z);
		/*
		 * Considering that z magnitude is the most likely to be the smaller, hence
		 * ensuring z <= y <= x, and not x <= y <= z, for less swaps.
		 */
		// Ensuring z <= y.
		if (z > y) {
			// y < z: swapping y and z
			double a = z;
			z = y;
			y = a;
		} else if (!(z <= y)) { // Testing if y or z is NaN.
			return hypot_NaN(x, y, z);
		}
		// Ensuring y <= x.
		if (z > x) {
			// x < z <= y: moving x
			double oldZ = z;
			z = x;
			double oldY = y;
			y = oldZ;
			x = oldY;
		} else if (y > x) {
			// z <= x < y: swapping x and y
			double a = y;
			y = x;
			x = a;
		} else if (x != x) { // Testing if x is NaN.
			return hypot_NaN(x, y, z);
		}

		if (x - y == x) {
			// y, hence z, too small to subtract from x.
			return x;
		} else if (y - z == y) {
			// z too small to subtract from y, hence x.
			double factor;
			if (x > HYPOT_MAX_MAG) {
				// x is too large: scaling down.
				x *= (1 / HYPOT_FACTOR);
				y *= (1 / HYPOT_FACTOR);
				factor = HYPOT_FACTOR;
			} else if (y < (1 / HYPOT_MAX_MAG)) {
				// y is too small: scaling up.
				x *= HYPOT_FACTOR;
				y *= HYPOT_FACTOR;
				factor = (1 / HYPOT_FACTOR);
			} else {
				factor = 1.0;
			}
			return factor * sqrt(x * x + y * y);
		} else {
			double factor;
			if (x > HYPOT_MAX_MAG) {
				// x is too large: scaling down.
				x *= (1 / HYPOT_FACTOR);
				y *= (1 / HYPOT_FACTOR);
				z *= (1 / HYPOT_FACTOR);
				factor = HYPOT_FACTOR;
			} else if (z < (1 / HYPOT_MAX_MAG)) {
				// z is too small: scaling up.
				x *= HYPOT_FACTOR;
				y *= HYPOT_FACTOR;
				z *= HYPOT_FACTOR;
				factor = (1 / HYPOT_FACTOR);
			} else {
				factor = 1.0;
			}
			// Adding smaller magnitudes together first.
			return factor * sqrt(x * x + (y * y + z * z));
		}
	}

	/*
	 * close values
	 */

	/**
	 * @param value A float value.
	 * @return Floor of value.
	 */
	public static float floor(float value) {
		final int exponent = getExponent(value);
		if (exponent < 0) {
			// abs(value) < 1.
			if (value < 0.0f) {
				return -1.0f;
			} else {
				// 0.0f, or -0.0f if value is -0.0f
				return 0.0f * value;
			}
		} else if (exponent < 23) {
			// A bit faster than using casts.
			final int bits = Float.floatToRawIntBits(value);
			final int anteCommaBits = bits & (0xFF800000 >> exponent);
			if ((value < 0.0f) && (anteCommaBits != bits)) {
				return Float.intBitsToFloat(anteCommaBits) - 1.0f;
			} else {
				return Float.intBitsToFloat(anteCommaBits);
			}
		} else {
			// +-Infinity, NaN, or a mathematical integer.
			return value;
		}
	}

	/**
	 * @param value A double value.
	 * @return Floor of value.
	 */
	public static double floor(double value) {
		if (USE_JDK_MATH) {
			return Math.floor(value);
		}
		if (ANTI_SLOW_CASTS) {
			double valueAbs = Math.abs(value);
			if (valueAbs <= (double) Integer.MAX_VALUE) {
				if (value > 0.0) {
					return (double) (int) value;
				} else if (value < 0.0) {
					double anteCommaDigits = (double) (int) value;
					if (value != anteCommaDigits) {
						return anteCommaDigits - 1.0;
					} else {
						return anteCommaDigits;
					}
				} else { // value is +-0.0 (not NaN due to test against Integer.MAX_VALUE)
					return value;
				}
			} else if (valueAbs < TWO_POW_52) {
				// We split the value in two:
				// high part, which is a mathematical integer,
				// and the rest, for which we can get rid of the
				// post comma digits by casting into an int.
				double highPart = ((int) (value * TWO_POW_N26)) * TWO_POW_26;
				if (value > 0.0) {
					return highPart + (double) ((int) (value - highPart));
				} else {
					double anteCommaDigits = highPart + (double) ((int) (value - highPart));
					if (value != anteCommaDigits) {
						return anteCommaDigits - 1.0;
					} else {
						return anteCommaDigits;
					}
				}
			} else { // abs(value) >= 2^52, or value is NaN
				return value;
			}
		} else {
			final int exponent = getExponent(value);
			if (exponent < 0) {
				// abs(value) < 1.
				if (value < 0.0) {
					return -1.0;
				} else {
					// 0.0, or -0.0 if value is -0.0
					return 0.0 * value;
				}
			} else if (exponent < 52) {
				// A bit faster than working on bits.
				final long matIntPart = (long) value;
				final double matIntToValue = value - (double) matIntPart;
				if (matIntToValue >= 0.0) {
					return (double) matIntPart;
				} else {
					return (double) (matIntPart - 1);
				}
			} else {
				// +-Infinity, NaN, or a mathematical integer.
				return value;
			}
		}
	}

	/**
	 * @param value A float value.
	 * @return Ceiling of value.
	 */
	public static float ceil(float value) {
		return -floor(-value);
	}

	/**
	 * @param value A double value.
	 * @return Ceiling of value.
	 */
	public static double ceil(double value) {
		if (USE_JDK_MATH) {
			return Math.ceil(value);
		}
		return -floor(-value);
	}

	/**
	 * Might have different semantics than Math.round(float), see bugs 6430675 and
	 * 8010430.
	 * 
	 * @param value A double value.
	 * @return Value rounded to nearest int, choosing superior int in case two are
	 *         equally close (i.e. rounding-up).
	 */
	public static int round(float value) {
		/*
		 * Not delegating to JDK, because we want delegation to provide at least as good
		 * results, and some supported JDK versions have bugged round() methods.
		 */
		// Algorithm by Dmitry Nadezhin (but replaced an if by a multiply)
		// (http://mail.openjdk.java.net/pipermail/core-libs-dev/2013-August/020247.html).
		final int bits = Float.floatToRawIntBits(value);
		final int biasedExp = ((bits >> 23) & 0xFF);
		// Shift to get rid of bits past comma except first one: will need to
		// 1-shift to the right to end up with correct magnitude.
		final int shift = (23 - 1 + MAX_FLOAT_EXPONENT) - biasedExp;
		if ((shift & -32) == 0) {
			int bitsSignum = (((bits >> 31) << 1) + 1);
			// shift in [0,31], so unbiased exp in [-9,22].
			int extendedMantissa = (0x00800000 | (bits & 0x007FFFFF)) * bitsSignum;
			// If value is positive and first bit past comma is 0, rounding
			// to lower integer, else to upper one, which is what "+1" and
			// then ">>1" do.
			return ((extendedMantissa >> shift) + 1) >> 1;
		} else {
			// +-Infinity, NaN, or a mathematical integer, or tiny.
			if (false && ANTI_SLOW_CASTS) { // not worth it
				if (Math.abs(value) >= -(float) Integer.MIN_VALUE) {
					// +-Infinity or a mathematical integer (mostly) out of int range.
					return (value < 0.0) ? Integer.MIN_VALUE : Integer.MAX_VALUE;
				}
				// NaN or a mathematical integer (mostly) in int range.
			}
			return (int) value;
		}
	}

	/**
	 * Might have different semantics than Math.round(double), see bugs 6430675 and
	 * 8010430.
	 * 
	 * @param value A double value.
	 * @return Value rounded to nearest long, choosing superior long in case two are
	 *         equally close (i.e. rounding-up).
	 */
	public static long round(double value) {
		/*
		 * Not delegating to JDK, because we want delegation to provide at least as good
		 * results, and some supported JDK versions have bugged round() methods.
		 */
		final long bits = Double.doubleToRawLongBits(value);
		final int biasedExp = (((int) (bits >> 52)) & 0x7FF);
		// Shift to get rid of bits past comma except first one: will need to
		// 1-shift to the right to end up with correct magnitude.
		final int shift = (52 - 1 + MAX_DOUBLE_EXPONENT) - biasedExp;
		if ((shift & -64) == 0) {
			long bitsSignum = (((bits >> 63) << 1) + 1);
			// shift in [0,63], so unbiased exp in [-12,51].
			long extendedMantissa = (0x0010000000000000L | (bits & 0x000FFFFFFFFFFFFFL)) * bitsSignum;
			// If value is positive and first bit past comma is 0, rounding
			// to lower integer, else to upper one, which is what "+1" and
			// then ">>1" do.
			return ((extendedMantissa >> shift) + 1L) >> 1;
		} else {
			// +-Infinity, NaN, or a mathematical integer, or tiny.
			if (ANTI_SLOW_CASTS) {
				if (Math.abs(value) >= -(double) Long.MIN_VALUE) {
					// +-Infinity or a mathematical integer (mostly) out of long range.
					return (value < 0.0) ? Long.MIN_VALUE : Long.MAX_VALUE;
				}
				// NaN or a mathematical integer (mostly) in long range.
			}
			return (long) value;
		}
	}

	/**
	 * @param value A float value.
	 * @return Value rounded to nearest int, choosing even int in case two are
	 *         equally close.
	 */
	public static int roundEven(float value) {
		final int sign = signFromBit(value);
		value = Math.abs(value);
		if (ANTI_SLOW_CASTS) {
			if (value < TWO_POW_23_F) {
				// Getting rid of post-comma bits.
				value = ((value + TWO_POW_23_F) - TWO_POW_23_F);
				return sign * (int) value;
			} else if (value < (float) Integer.MAX_VALUE) { // "<=" doesn't work, because of float precision
				// value is in [-Integer.MAX_VALUE,Integer.MAX_VALUE]
				return sign * (int) value;
			}
		} else {
			if (value < TWO_POW_23_F) {
				// Getting rid of post-comma bits.
				value = ((value + TWO_POW_23_F) - TWO_POW_23_F);
			}
		}
		return (int) (sign * value);
	}

	/**
	 * @param value A double value.
	 * @return Value rounded to nearest long, choosing even long in case two are
	 *         equally close.
	 */
	public static long roundEven(double value) {
		final int sign = (int) signFromBit(value);
		value = Math.abs(value);
		if (value < TWO_POW_52) {
			// Getting rid of post-comma bits.
			value = ((value + TWO_POW_52) - TWO_POW_52);
		}
		if (ANTI_SLOW_CASTS) {
			if (value <= (double) Integer.MAX_VALUE) {
				// value is in [-Integer.MAX_VALUE,Integer.MAX_VALUE]
				return sign * (int) value;
			}
		}
		return (long) (sign * value);
	}

	/**
	 * @param value A float value.
	 * @return The float mathematical integer closest to the specified value,
	 *         choosing even one if two are equally close, or respectively NaN,
	 *         +-Infinity or +-0.0f if the value is any of these.
	 */
	public static float rint(float value) {
		final int sign = signFromBit(value);
		value = Math.abs(value);
		if (value < TWO_POW_23_F) {
			// Getting rid of post-comma bits.
			value = ((TWO_POW_23_F + value) - TWO_POW_23_F);
		}
		// Restoring original sign.
		return sign * value;
	}

	/**
	 * @param value A double value.
	 * @return The double mathematical integer closest to the specified value,
	 *         choosing even one if two are equally close, or respectively NaN,
	 *         +-Infinity or +-0.0 if the value is any of these.
	 */
	public static double rint(double value) {
		if (USE_JDK_MATH) {
			return Math.rint(value);
		}
		final int sign = (int) signFromBit(value);
		value = Math.abs(value);
		if (value < TWO_POW_52) {
			// Getting rid of post-comma bits.
			value = ((TWO_POW_52 + value) - TWO_POW_52);
		}
		// Restoring original sign.
		return sign * value;
	}

	/*
	 * close int values
	 * 
	 * Never delegating to JDK for these methods, for we should always be faster and
	 * exact, and JDK doesn't exactly have such methods.
	 */

	/**
	 * @param value A double value.
	 * @return Floor of value as int, or closest int if floor is out of int range,
	 *         or 0 if value is NaN.
	 */
	public static int floorToInt(double value) {
		int valueInt = (int) value;
		if (value < 0.0) {
			if (value == (double) valueInt) {
				return valueInt;
			} else {
				if (valueInt == Integer.MIN_VALUE) {
					return valueInt;
				} else {
					return valueInt - 1;
				}
			}
		} else { // >= 0 or NaN.
			return valueInt;
		}
	}

	/**
	 * @param value A double value.
	 * @return Ceiling of value as int, or closest int if ceiling is out of int
	 *         range, or 0 if value is NaN.
	 */
	public static int ceilToInt(double value) {
		int valueInt = (int) value;
		if (value > 0.0) {
			if (value == (double) valueInt) {
				return valueInt;
			} else {
				if (valueInt == Integer.MAX_VALUE) {
					return valueInt;
				} else {
					return valueInt + 1;
				}
			}
		} else { // <= 0 or NaN.
			return valueInt;
		}
	}

	/**
	 * @param value A double value.
	 * @return Value rounded to nearest int, choosing superior int in case two are
	 *         equally close (i.e. rounding-up).
	 */
	public static int roundToInt(double value) {
		/*
		 * We don't gain much by reimplementing rounding, except for pathologically
		 * large values, which should not be a common case when dealing with ints, so we
		 * just use round(double).
		 */
		return NumbersUtils.toInt(round(value));
	}

	/**
	 * @param value A double value.
	 * @return Value rounded to nearest int, choosing even int in case two are
	 *         equally close.
	 */
	public static int roundEvenToInt(double value) {
		final int sign = (int) signFromBit(value);
		value = Math.abs(value);
		/*
		 * Applying the post-comma bits removal logic even if value is out of int range,
		 * to avoid a test, for it doesn't mess up the result, and we want to optimize
		 * for the case of values in int range.
		 */
		value = ((value + TWO_POW_52) - TWO_POW_52);
		return (int) (sign * value);
	}

	/*
	 * ranges
	 */

	/**
	 * @param min   A float value.
	 * @param max   A float value.
	 * @param value A float value.
	 * @return min if value < min, max if value > max, value otherwise.
	 */
	public static float toRange(float min, float max, float value) {
		return NumbersUtils.toRange(min, max, value);
	}

	/**
	 * @param min   A double value.
	 * @param max   A double value.
	 * @param value A double value.
	 * @return min if value < min, max if value > max, value otherwise.
	 */
	public static double toRange(double min, double max, double value) {
		return NumbersUtils.toRange(min, max, value);
	}

	/*
	 * binary operators (/,%)
	 */

	/**
	 * Returns dividend - divisor * n, where n is the mathematical integer closest
	 * to dividend/divisor. If dividend/divisor is equally close to surrounding
	 * integers, we choose n to be the integer of smallest magnitude, which makes
	 * this treatment differ from Math.IEEEremainder(double,double), where n is
	 * chosen to be the even integer. Note that the choice of n is not done
	 * considering the double approximation of dividend/divisor, because it could
	 * cause result to be outside [-|divisor|/2,|divisor|/2] range. The practical
	 * effect is that if multiple results would be possible, we always choose the
	 * result that is the closest to (and has the same sign as) the dividend. Ex. :
	 * - for (-3.0,2.0), this method returns -1.0, whereas Math.IEEEremainder
	 * returns 1.0. - for (-5.0,2.0), both this method and Math.IEEEremainder return
	 * -1.0.
	 * 
	 * If the remainder is zero, its sign is the same as the sign of the first
	 * argument. If either argument is NaN, or the first argument is infinite, or
	 * the second argument is positive zero or negative zero, then the result is
	 * NaN. If the first argument is finite and the second argument is infinite,
	 * then the result is the same as the first argument.
	 * 
	 * NB: - Modulo operator (%) returns a value in ]-|divisor|,|divisor|[, which
	 * sign is the same as dividend. - As for modulo operator, the sign of the
	 * divisor has no effect on the result. - On some architecture, % operator has
	 * been observed to return NaN for some subnormal values of divisor, when
	 * dividend exponent is 1023, which impacts the correctness of this method.
	 * 
	 * @param dividend Dividend.
	 * @param divisor  Divisor.
	 * @return Remainder of dividend/divisor, i.e. a value in
	 *         [-|divisor|/2,|divisor|/2].
	 */
	public static double remainder(double dividend, double divisor) {
		if (Double.isInfinite(divisor)) {
			if (Double.isInfinite(dividend)) {
				return Double.NaN;
			} else {
				return dividend;
			}
		}
		double value = dividend % divisor;
		if (Math.abs(value + value) > Math.abs(divisor)) {
			return value + ((value > 0.0) ? -Math.abs(divisor) : Math.abs(divisor));
		} else {
			return value;
		}
	}

	/**
	 * @param angle Angle in radians.
	 * @return The same angle, in radians, but in [-PI,PI].
	 */
	public static double normalizeMinusPiPi(double angle) {
		// Not modifying values in output range.
		if ((angle >= -Math.PI) && (angle <= Math.PI)) {
			return angle;
		}
		return remainderTwoPi(angle);
	}

	/**
	 * Not accurate for large values.
	 * 
	 * @param angle Angle in radians.
	 * @return The same angle, in radians, but in [-PI,PI].
	 */
	public static double normalizeMinusPiPiFast(double angle) {
		// Not modifying values in output range.
		if ((angle >= -Math.PI) && (angle <= Math.PI)) {
			return angle;
		}
		return remainderTwoPiFast(angle);
	}

	/**
	 * @param angle Angle in radians.
	 * @return The same angle, in radians, but in [0,2*PI].
	 */
	public static double normalizeZeroTwoPi(double angle) {
		// Not modifying values in output range.
		if ((angle >= 0.0) && (angle <= 2 * Math.PI)) {
			return angle;
		}
		angle = remainderTwoPi(angle);
		if (angle < 0.0) {
			// LO then HI is theoretically better (when starting near 0).
			return (angle + TWOPI_LO) + TWOPI_HI;
		} else {
			return angle;
		}
	}

	/**
	 * Not accurate for large values.
	 * 
	 * @param angle Angle in radians.
	 * @return The same angle, in radians, but in [0,2*PI].
	 */
	public static double normalizeZeroTwoPiFast(double angle) {
		// Not modifying values in output range.
		if ((angle >= 0.0) && (angle <= 2 * Math.PI)) {
			return angle;
		}
		angle = remainderTwoPiFast(angle);
		if (angle < 0.0) {
			// LO then HI is theoretically better (when starting near 0).
			return (angle + TWOPI_LO) + TWOPI_HI;
		} else {
			return angle;
		}
	}

	/**
	 * @param angle Angle in radians.
	 * @return Angle value modulo PI, in radians, in [-PI/2,PI/2].
	 */
	public static double normalizeMinusHalfPiHalfPi(double angle) {
		// Not modifying values in output range.
		if ((angle >= -Math.PI / 2) && (angle <= Math.PI / 2)) {
			return angle;
		}
		return remainderPi(angle);
	}

	/**
	 * Not accurate for large values.
	 * 
	 * @param angle Angle in radians.
	 * @return Angle value modulo PI, in radians, in [-PI/2,PI/2].
	 */
	public static double normalizeMinusHalfPiHalfPiFast(double angle) {
		// Not modifying values in output range.
		if ((angle >= -Math.PI / 2) && (angle <= Math.PI / 2)) {
			return angle;
		}
		return remainderPiFast(angle);
	}

	/*
	 * floating points utils
	 */

	/**
	 * @param value A float value.
	 * @return true if the specified value is NaN or +-Infinity, false otherwise.
	 */
	public static boolean isNaNOrInfinite(float value) {
		return NumbersUtils.isNaNOrInfinite(value);
	}

	/**
	 * @param value A double value.
	 * @return true if the specified value is NaN or +-Infinity, false otherwise.
	 */
	public static boolean isNaNOrInfinite(double value) {
		return NumbersUtils.isNaNOrInfinite(value);
	}

	/**
	 * @param value A float value.
	 * @return Value unbiased exponent.
	 */
	public static int getExponent(float value) {
		return ((Float.floatToRawIntBits(value) >> 23) & 0xFF) - MAX_FLOAT_EXPONENT;
	}

	/**
	 * @param value A double value.
	 * @return Value unbiased exponent.
	 */
	public static int getExponent(double value) {
		return (((int) (Double.doubleToRawLongBits(value) >> 52)) & 0x7FF) - MAX_DOUBLE_EXPONENT;
	}

	/**
	 * @param value A float value.
	 * @return -1.0f if the specified value is < 0, 1.0f if it is > 0, and the value
	 *         itself if it is NaN or +-0.0f.
	 */
	public static float signum(float value) {
		if (USE_JDK_MATH) {
			return Math.signum(value);
		}
		if ((value == 0.0f) || (value != value)) {
			return value;
		}
		return (float) signFromBit(value);
	}

	/**
	 * @param value A double value.
	 * @return -1.0 if the specified value is < 0, 1.0 if it is > 0, and the value
	 *         itself if it is NaN or +-0.0.
	 */
	public static double signum(double value) {
		if (USE_JDK_MATH) {
			return Math.signum(value);
		}
		if ((value == 0.0) || (value != value)) {
			return value;
		}
		if (ANTI_SLOW_CASTS) {
			return (double) (int) signFromBit(value);
		} else {
			return (double) signFromBit(value);
		}
	}

	/**
	 * @param value A float value.
	 * @return -1 if sign bit is 1, 1 if sign bit is 0.
	 */
	public static int signFromBit(float value) {
		return ((Float.floatToRawIntBits(value) >> 30) | 1);
	}

	/**
	 * @param value A double value.
	 * @return -1 if sign bit is 1, 1 if sign bit is 0.
	 */
	public static long signFromBit(double value) {
		// Returning a long, to avoid useless cast into int.
		return ((Double.doubleToRawLongBits(value) >> 62) | 1);
	}

	/**
	 * A sign of NaN can be interpreted as positive or negative.
	 *
	 * @param magnitude A float value.
	 * @param sign      A float value.
	 * @return A value with the magnitude of the first argument, and the sign of the
	 *         second argument.
	 */
	public static float copySign(float magnitude, float sign) {
		return Float.intBitsToFloat((Float.floatToRawIntBits(sign) & Integer.MIN_VALUE) | (Float.floatToRawIntBits(magnitude) & Integer.MAX_VALUE));
	}

	/**
	 * A sign of NaN can be interpreted as positive or negative.
	 *
	 * @param magnitude A double value.
	 * @param sign      A double value.
	 * @return A value with the magnitude of the first argument, and the sign of the
	 *         second argument.
	 */
	public static double copySign(double magnitude, double sign) {
		return Double.longBitsToDouble((Double.doubleToRawLongBits(sign) & Long.MIN_VALUE) | (Double.doubleToRawLongBits(magnitude) & Long.MAX_VALUE));
	}

	/**
	 * The ULP (Unit in the Last Place) is the distance to the next value larger in
	 * magnitude.
	 *
	 * @param value A float value.
	 * @return The size of an ulp of the specified value, or Float.MIN_VALUE if it
	 *         is +-0.0f, or +Infinity if it is +-Infinity, or NaN if it is NaN.
	 */
	public static float ulp(float value) {
		if (USE_JDK_MATH) {
			return Math.ulp(value);
		}
		/*
		 * Look-up table not really worth it in micro-benchmark, so should be worse with
		 * cache-misses.
		 */
		final int exponent = getExponent(value);
		if (exponent >= (MIN_FLOAT_NORMAL_EXPONENT + 23)) {
			if (exponent == MAX_FLOAT_EXPONENT + 1) {
				// NaN or +-Infinity
				return Math.abs(value);
			}
			// normal: returning 2^(exponent-23)
			return Float.intBitsToFloat((exponent + (MAX_FLOAT_EXPONENT - 23)) << 23);
		} else {
			if (exponent == MIN_FLOAT_NORMAL_EXPONENT - 1) {
				// +-0.0f or subnormal
				return Float.MIN_VALUE;
			}
			// subnormal result
			return Float.intBitsToFloat(1 << (exponent - MIN_FLOAT_NORMAL_EXPONENT));
		}
	}

	/**
	 * The ULP (Unit in the Last Place) is the distance to the next value larger in
	 * magnitude.
	 *
	 * @param value A double value.
	 * @return The size of an ulp of the specified value, or Double.MIN_VALUE if it
	 *         is +-0.0, or +Infinity if it is +-Infinity, or NaN if it is NaN.
	 */
	public static double ulp(double value) {
		if (USE_JDK_MATH) {
			return Math.ulp(value);
		}
		/*
		 * Look-up table not really worth it in micro-benchmark, so should be worse with
		 * cache-misses.
		 */
		final int exponent = getExponent(value);
		if (exponent >= (MIN_DOUBLE_NORMAL_EXPONENT + 52)) {
			if (exponent == MAX_DOUBLE_EXPONENT + 1) {
				// NaN or +-Infinity
				return Math.abs(value);
			}
			// normal: returning 2^(exponent-52)
			return Double.longBitsToDouble((exponent + (MAX_DOUBLE_EXPONENT - 52L)) << 52);
		} else {
			if (exponent == MIN_DOUBLE_NORMAL_EXPONENT - 1) {
				// +-0.0f or subnormal
				return Double.MIN_VALUE;
			}
			// subnormal result
			return Double.longBitsToDouble(1L << (exponent - MIN_DOUBLE_NORMAL_EXPONENT));
		}
	}

	/**
	 * If both arguments are +-0.0(f), (float)direction is returned.
	 * 
	 * If both arguments are +Infinity or -Infinity, respectively +Infinity or
	 * -Infinity is returned.
	 *
	 * @param start     A float value.
	 * @param direction A double value.
	 * @return The float adjacent to start towards direction, considering that
	 *         +(-)Float.MIN_VALUE is adjacent to +(-)0.0f, and that
	 *         +(-)Float.MAX_VALUE is adjacent to +(-)Infinity, or NaN if any
	 *         argument is NaN.
	 */
	public static float nextAfter(float start, double direction) {
		if (direction < start) {
			// Going towards -Infinity.
			if (start == 0.0f) {
				// +-0.0f
				return -Float.MIN_VALUE;
			}
			final int bits = Float.floatToRawIntBits(start);
			return Float.intBitsToFloat(bits + ((bits > 0) ? -1 : 1));
		} else if (direction > start) {
			// Going towards +Infinity.
			// +0.0f to get rid of eventual -0.0f
			final int bits = Float.floatToRawIntBits(start + 0.0f);
			return Float.intBitsToFloat(bits + (bits >= 0 ? 1 : -1));
		} else if (start == direction) {
			return (float) direction;
		} else {
			// Returning a NaN derived from the input NaN(s).
			return start + (float) direction;
		}
	}

	/**
	 * If both arguments are +-0.0, direction is returned.
	 * 
	 * If both arguments are +Infinity or -Infinity, respectively +Infinity or
	 * -Infinity is returned.
	 *
	 * @param start     A double value.
	 * @param direction A double value.
	 * @return The double adjacent to start towards direction, considering that
	 *         +(-)Double.MIN_VALUE is adjacent to +(-)0.0, and that
	 *         +(-)Double.MAX_VALUE is adjacent to +(-)Infinity, or NaN if any
	 *         argument is NaN.
	 */
	public static double nextAfter(double start, double direction) {
		if (direction < start) {
			// Going towards -Infinity.
			if (start == 0.0) {
				// +-0.0
				return -Double.MIN_VALUE;
			}
			final long bits = Double.doubleToRawLongBits(start);
			return Double.longBitsToDouble(bits + ((bits > 0) ? -1 : 1));
		} else if (direction > start) {
			// Going towards +Infinity.
			// +0.0 to get rid of eventual -0.0
			final long bits = Double.doubleToRawLongBits(start + 0.0f);
			return Double.longBitsToDouble(bits + (bits >= 0 ? 1 : -1));
		} else if (start == direction) {
			return direction;
		} else {
			// Returning a NaN derived from the input NaN(s).
			return start + direction;
		}
	}

	/**
	 * Semantically equivalent to nextAfter(start,Double.NEGATIVE_INFINITY).
	 */
	public static float nextDown(float start) {
		if (start > Float.NEGATIVE_INFINITY) {
			if (start == 0.0f) {
				// +-0.0f
				return -Float.MIN_VALUE;
			}
			final int bits = Float.floatToRawIntBits(start);
			return Float.intBitsToFloat(bits + ((bits > 0) ? -1 : 1));
		} else if (start == Float.NEGATIVE_INFINITY) {
			return Float.NEGATIVE_INFINITY;
		} else {
			// NaN
			return start;
		}
	}

	/**
	 * Semantically equivalent to nextAfter(start,Double.NEGATIVE_INFINITY).
	 */
	public static double nextDown(double start) {
		if (start > Double.NEGATIVE_INFINITY) {
			if (start == 0.0) {
				// +-0.0
				return -Double.MIN_VALUE;
			}
			final long bits = Double.doubleToRawLongBits(start);
			return Double.longBitsToDouble(bits + ((bits > 0) ? -1 : 1));
		} else if (start == Double.NEGATIVE_INFINITY) {
			return Double.NEGATIVE_INFINITY;
		} else {
			// NaN
			return start;
		}
	}

	/**
	 * Semantically equivalent to nextAfter(start,Double.POSITIVE_INFINITY).
	 */
	public static float nextUp(float start) {
		if (start < Float.POSITIVE_INFINITY) {
			// +0.0f to get rid of eventual -0.0f
			final int bits = Float.floatToRawIntBits(start + 0.0f);
			return Float.intBitsToFloat(bits + (bits >= 0 ? 1 : -1));
		} else if (start == Float.POSITIVE_INFINITY) {
			return Float.POSITIVE_INFINITY;
		} else {
			// NaN
			return start;
		}
	}

	/**
	 * Semantically equivalent to nextAfter(start,Double.POSITIVE_INFINITY).
	 */
	public static double nextUp(double start) {
		if (start < Double.POSITIVE_INFINITY) {
			// +0.0 to get rid of eventual -0.0
			final long bits = Double.doubleToRawLongBits(start + 0.0);
			return Double.longBitsToDouble(bits + (bits >= 0 ? 1 : -1));
		} else if (start == Double.POSITIVE_INFINITY) {
			return Double.POSITIVE_INFINITY;
		} else {
			// NaN
			return start;
		}
	}

	/**
	 * Precision may be lost if the result is subnormal.
	 *
	 * @param value       A float value.
	 * @param scaleFactor An int value.
	 * @return value * 2^scaleFactor, or a value equivalent to the specified one if
	 *         it is NaN, +-Infinity or +-0.0f.
	 */
	public static float scalb(float value, int scaleFactor) {
		// Large enough to imply overflow or underflow for
		// a finite non-zero value.
		final int MAX_SCALE = 2 * MAX_FLOAT_EXPONENT + 23 + 1;

		// Making sure scaling factor is in a reasonable range.
		scaleFactor = Math.max(Math.min(scaleFactor, MAX_SCALE), -MAX_SCALE);

		return (float) (((double) value) * twoPowNormal(scaleFactor));
	}

	/**
	 * Precision may be lost if the result is subnormal.
	 *
	 * @param value       A double value.
	 * @param scaleFactor An int value.
	 * @return value * 2^scaleFactor, or a value equivalent to the specified one if
	 *         it is NaN, +-Infinity or +-0.0.
	 */
	public static double scalb(double value, int scaleFactor) {
		if ((scaleFactor > -MAX_DOUBLE_EXPONENT) && (scaleFactor <= MAX_DOUBLE_EXPONENT)) {
			// Quick case (as done in apache FastMath).
			return value * twoPowNormal(scaleFactor);
		}

		// Large enough to imply overflow or underflow for
		// a finite non-zero value.
		final int MAX_SCALE = 2 * MAX_DOUBLE_EXPONENT + 52 + 1;

		// Making sure scaling factor is in a reasonable range.
		final int exponentAdjust;
		final int scaleIncrement;
		final double exponentDelta;
		if (scaleFactor < 0) {
			scaleFactor = Math.max(scaleFactor, -MAX_SCALE);
			scaleIncrement = -512;
			exponentDelta = TWO_POW_N512;
		} else {
			scaleFactor = Math.min(scaleFactor, MAX_SCALE);
			scaleIncrement = 512;
			exponentDelta = TWO_POW_512;
		}

		// Calculating (scaleFactor % +-512), 512 = 2^9, using
		// technique from "Hacker's Delight" section 10-2.
		final int t = ((scaleFactor >> (9 - 1)) >>> (32 - 9));
		exponentAdjust = ((scaleFactor + t) & (512 - 1)) - t;

		value *= twoPowNormal(exponentAdjust);
		scaleFactor -= exponentAdjust;

		while (scaleFactor != 0) {
			value *= exponentDelta;
			scaleFactor -= scaleIncrement;
		}

		return value;
	}

	/*
	 * Non-redefined Math public values and treatments.
	 */

	public static float abs(float a) {
		return Math.abs(a);
	}

	public static double abs(double a) {
		return Math.abs(a);
	}

	public static float min(float a, float b) {
		return Math.min(a, b);
	}

	public static double min(double a, double b) {
		return Math.min(a, b);
	}

	public static float max(float a, float b) {
		return Math.max(a, b);
	}

	public static double max(double a, double b) {
		return Math.max(a, b);
	}

	public static double IEEEremainder(double f1, double f2) {
		return Math.IEEEremainder(f1, f2);
	}

	public static double random() {
		return RandomUtil.random();
	}

	// --------------------------------------------------------------------------
	// PRIVATE METHODS
	// --------------------------------------------------------------------------

	/**
	 * Instantiable.
	 */
	public Jafama() {
	}

	/*
	 * Remainders (accurate).
	 */

	/**
	 * @param angle Angle in radians.
	 * @return Remainder of (angle % (2*PI)), in [-PI,PI].
	 */
	private static double remainderTwoPi(double angle) {
		if (USE_JDK_MATH) {
			return jdkRemainderTwoPi(angle);
		}
		boolean negateResult = false;
		if (angle < 0.0) {
			angle = -angle;
			negateResult = true;
		}
		if (angle <= (4 * NORMALIZE_ANGLE_MAX_MEDIUM_DOUBLE_PIO2)) {
			double fn = (double) (int) (angle * TWOPI_INV + 0.5);
			angle = (angle - fn * TWOPI_HI) - fn * TWOPI_LO;
			// Ensuring range.
			// HI/LO can help a bit, even though we are always far from 0.
			if (angle < -Math.PI) {
				angle = (angle + TWOPI_HI) + TWOPI_LO;
			} else if (angle > Math.PI) {
				angle = (angle - TWOPI_HI) - TWOPI_LO;
			}
			return negateResult ? -angle : angle;
		} else if (angle < Double.POSITIVE_INFINITY) {
			angle = heavyRemainderTwoPi(angle);
			return negateResult ? -angle : angle;
		} else { // angle is +Infinity or NaN
			return Double.NaN;
		}
	}

	/**
	 * @param angle Angle in radians.
	 * @return Remainder of (angle % PI), in [-PI/2,PI/2].
	 */
	private static double remainderPi(double angle) {
		if (USE_JDK_MATH) {
			return jdkRemainderPi(angle);
		}
		boolean negateResult = false;
		if (angle < 0.0) {
			angle = -angle;
			negateResult = true;
		}
		if (angle <= (2 * NORMALIZE_ANGLE_MAX_MEDIUM_DOUBLE_PIO2)) {
			double fn = (double) (int) (angle * PI_INV + 0.5);
			angle = (angle - fn * PI_HI) - fn * PI_LO;
			// Ensuring range.
			// HI/LO can help a bit, even though we are always far from 0.
			if (angle < -Math.PI / 2) {
				angle = (angle + PI_HI) + PI_LO;
			} else if (angle > Math.PI / 2) {
				angle = (angle - PI_HI) - PI_LO;
			}
			return negateResult ? -angle : angle;
		} else if (angle < Double.POSITIVE_INFINITY) {
			angle = heavyRemainderPi(angle);
			return negateResult ? -angle : angle;
		} else { // angle is +Infinity or NaN
			return Double.NaN;
		}
	}

	/**
	 * @param angle Angle in radians.
	 * @return Bits of double corresponding to remainder of (angle % (PI/2)), in
	 *         [-PI/4,PI/4], with quadrant encoded in exponent bits.
	 */
	private static long remainderPiO2(double angle) {
		if (USE_JDK_MATH) {
			return jdkRemainderPiO2(angle, false);
		}
		boolean negateResult = false;
		if (angle < 0.0) {
			angle = -angle;
			negateResult = true;
		}
		if (angle <= NORMALIZE_ANGLE_MAX_MEDIUM_DOUBLE_PIO2) {
			int n = (int) (angle * PIO2_INV + 0.5);
			double fn = (double) n;
			angle = (angle - fn * PIO2_HI) - fn * PIO2_LO;
			// Ensuring range.
			// HI/LO can help a bit, even though we are always far from 0.
			if (angle < -Math.PI / 4) {
				angle = (angle + PIO2_HI) + PIO2_LO;
				n--;
			} else if (angle > Math.PI / 4) {
				angle = (angle - PIO2_HI) - PIO2_LO;
				n++;
			}
			if (negateResult) {
				angle = -angle;
			}
			return encodeRemainderAndQuadrant(angle, n & 3);
		} else if (angle < Double.POSITIVE_INFINITY) {
			return heavyRemainderPiO2(angle, negateResult);
		} else { // angle is +Infinity or NaN
			return encodeRemainderAndQuadrant(Double.NaN, 0);
		}
	}

	/*
	 * Remainders (fast).
	 */

	/**
	 * Not accurate for large values.
	 * 
	 * @param angle Angle in radians.
	 * @return Remainder of (angle % (2*PI)), in [-PI,PI].
	 */
	private static double remainderTwoPiFast(double angle) {
		if (USE_JDK_MATH) {
			return jdkRemainderTwoPi(angle);
		}
		boolean negateResult = false;
		if (angle < 0.0) {
			angle = -angle;
			negateResult = true;
		}
		// - We don't bother with values higher than (2*PI*(2^52)),
		// since they are spaced by 2*PI or more from each other.
		// - For large values, we don't use % because it might be very slow,
		// and we split computation in two, because cast from double to int
		// with large numbers might be very slow also.
		if (angle <= TWO_POW_26 * (2 * Math.PI)) {
			// ok
		} else if (angle <= TWO_POW_52 * (2 * Math.PI)) {
			// Computing remainder of angle modulo TWO_POW_26*(2*PI).
			double fn = (double) (int) (angle * (TWOPI_INV / TWO_POW_26) + 0.5);
			angle = (angle - fn * (TWOPI_HI * TWO_POW_26)) - fn * (TWOPI_LO * TWO_POW_26);
			// Here, angle is in [-TWO_POW_26*PI,TWO_POW_26*PI], or so.
			if (angle < 0.0) {
				angle = -angle;
				negateResult = !negateResult;
			}
		} else if (angle < Double.POSITIVE_INFINITY) {
			return 0.0;
		} else { // angle is +Infinity or NaN
			return Double.NaN;
		}

		// Computing remainder of angle modulo 2*PI.
		double fn = (double) (int) (angle * TWOPI_INV + 0.5);
		angle = (angle - fn * TWOPI_HI) - fn * TWOPI_LO;

		// Ensuring range.
		// HI/LO can help a bit, even though we are always far from 0.
		if (angle < -Math.PI) {
			angle = (angle + TWOPI_HI) + TWOPI_LO;
		} else if (angle > Math.PI) {
			angle = (angle - TWOPI_HI) - TWOPI_LO;
		}
		return negateResult ? -angle : angle;
	}

	/**
	 * Not accurate for large values.
	 * 
	 * @param angle Angle in radians.
	 * @return Remainder of (angle % PI), in [-PI/2,PI/2].
	 */
	private static double remainderPiFast(double angle) {
		if (USE_JDK_MATH) {
			return jdkRemainderPi(angle);
		}
		boolean negateResult = false;
		if (angle < 0.0) {
			angle = -angle;
			negateResult = true;
		}
		// - We don't bother with values higher than (PI*(2^52)),
		// since they are spaced by PI or more from each other.
		// - For large values, we don't use % because it might be very slow,
		// and we split computation in two, because cast from double to int
		// with large numbers might be very slow also.
		if (angle <= TWO_POW_26 * Math.PI) {
			// ok
		} else if (angle <= TWO_POW_52 * Math.PI) {
			// Computing remainder of angle modulo TWO_POW_26*PI.
			double fn = (double) (int) (angle * (PI_INV / TWO_POW_26) + 0.5);
			angle = (angle - fn * (PI_HI * TWO_POW_26)) - fn * (PI_LO * TWO_POW_26);
			// Here, angle is in [-TWO_POW_26*PI/2,TWO_POW_26*PI/2], or so.
			if (angle < 0.0) {
				angle = -angle;
				negateResult = !negateResult;
			}
		} else if (angle < Double.POSITIVE_INFINITY) {
			return 0.0;
		} else { // angle is +Infinity or NaN
			return Double.NaN;
		}

		// Computing remainder of angle modulo PI.
		double fn = (double) (int) (angle * PI_INV + 0.5);
		angle = (angle - fn * PI_HI) - fn * PI_LO;

		// Ensuring range.
		// HI/LO can help a bit, even though we are always far from 0.
		if (angle < -Math.PI / 2) {
			angle = (angle + PI_HI) + PI_LO;
		} else if (angle > Math.PI / 2) {
			angle = (angle - PI_HI) - PI_LO;
		}
		return negateResult ? -angle : angle;
	}
}